Challenges of Big Data Analysis

Big Data bring new opportunities to modern society and challenges to data scientists. On one hand, Big Data hold great promises for discovering subtle population patterns and heterogeneities that are not possible with small-scale data. On the other hand, the massive sample size and high dimensionality of Big Data introduce unique computational and statistical challenges, including scalability and storage bottleneck, noise accumulation, spurious correlation, incidental endogeneity, and measurement errors. These challenges are distinguished and require new computational and statistical paradigm. This article give overviews on the salient features of Big Data and how these features impact on paradigm change on statistical and computational methods as well as computing architectures. We also provide various new perspectives on the Big Data analysis and computation. In particular, we emphasis on the viability of the sparsest solution in high-confidence set and point out that exogeneous assumptions in most statistical methods for Big Data can not be validated due to incidental endogeneity. They can lead to wrong statistical inferences and consequently wrong scientific conclusions.


Introduction
Big Data promise new levels of scientific discovery and economic value.What is new about Big Data and how they differ from the traditional small or medium-scale data?This article overviews the opportunities and challenges brought by Big Data, with emphasis on the distinguished features of Big Data and statistical and computational methods as well as computing architecture to deal with them.

Background
We are entering the era of Big Data -a term that refers to the explosion of available information.Such a Big Data movement is driven by the fact that massive amounts of very high dimensional or unstructured data are continuously produced and stored with much cheaper cost than they used to be.For example, in genomics we have seen a dramatic drop in price for whole genome sequencing (Stein, 2010).This is also true in other areas like social media analysis, biomedical imaging, high frequency finance, analysis of surveillance videos and retail sales.The existing trend that data can be produced and stored more massively and cheaply is likely to maintain or even accelerate in the future (Donoho, 2000).This trend will have deep impact on science, engineering, and business.For example, scientific advances are becoming more and more data-driven and researchers will more and more think of themselves as consumers of data.The massive amounts of high dimensional data bring both opportunities and new challenges to data analysis.Valid statistical analysis for Big Data is becoming increasingly important.

Goals and Challenges of Analyzing Big Data
What are the goals of analyzing Big Data?According to Bickel (2008), two main goals of high-dimensional data analysis are to develop effective methods that can accurately predict the future observations and at the same time to gain insight into the relationship between the features and response for scientific purposes.Furthermore, due to large sample size, Big Data give rise to two additional goals: To understand heterogeneity and commonality across different sub-populations.In other word, Big Data give promises for: (i) Exploring the hidden structures of each sub-population of the data, which is traditionally not feasible and might even be treated as "outliers" when the sample size is small; (ii) Extracting important common features across many sub-populations even when there are large individual variations.
What are the challenges of analyzing Big Data?Big Data are characterized by high dimensionality and large sample size.These two features raise three unique challenges: (i) High dimensionality brings noise accumulation, spurious correlations, and incidental homogeneity; (ii) High dimensionality combined with large sample size creates issues such as heavy computational cost and algorithmic instability; (iii) The massive samples in Big 1998; Wu and Lange, 2008), iterative shrinkage-thresholding algorithms (Daubechies et al., 2004;Beck and Teboulle, 2009), are proposed.Besides large-scale optimization algorithms, Big Data also motivate the development of majorization-minimization algorithms (Lange et al., 2000;Hunter and Li, 2005;Zou and Li, 2008),"large-scale screening and small-scale optimization" framework (Fan et al., 2009), parallel computing methods (Boyd et al., 2011;Bradley et al., 2011;Low et al., 2012), and approximate algorithms that are scalable to large sample size.

Organization of This Article
The rest of this article is organized as follows.Section 2 overviews the rise of Big Data problem from science, engineering, and social science.Section 3 explains some unique features of Big Data and their impacts on statistical inference.Statistical methods that tackle these Big Data problems are given in Section 4. Section 5 gives an overview on scalable computing infrastructure for Big Data storage and processing.Section 6 discusses the computational aspect of Big Data and introduces some recent progress.Section 7 concludes.

Rises of Big Data
Massive sample size and high dimensionality characterize many contemporary datasets.For example, in genomics, there have been more than 500, 000 microarrays that are publicly available with each array containing tens of thousands of expression values of molecules; In biomedical engineering, there have been tens of thousands of terabytes of fMRI images with each image containing more than 50, 000 voxel values.Other examples of massive and high dimensional data include unstructured text corpus, social medias, financial time series, e-commerce data, retail transaction records, surveillance videos.We now briefly illustrate some of these Big Data problems.

Genomics
Many new technologies have been developed in genomics and enable inexpensive and highthroughput measurement of the whole genome and transcriptome.These technologies allow biologists to generate hundreds of thousands of datasets and have shifted their primary interests from the acquisition of biological sequences to the study of biological function.The availability of massive datasets sheds light towards new scientific discoveries.For example, the large amount of genome sequencing data now make it possible to uncover the genetic markers of rare disorders (Worthey et al., 2010;Chen et al., 2012) and find associations between diseases and rare sequence variants (Cohen et al., 2004;Han and Pan, 2010).The breakthroughs in biomedical imaging technology allow scientists to simultaneously monitor many gene and protein functions, permitting us to study interactions in regulatory processes and neuron activities.Moreover, the emergence of publicly available genomic databases enables integrative analysis which combines information from many sources for drawing scientific conclusions.These researches give rise to many computational methods as well as new statistical thinking and challenges (Bickel et al., 2009a).
One of the important steps in genomic data analysis is to remove systematic biases (e.g., intensity effect, batch effect, dye effect, block effect, among others).Such systematic biases are due to experimental variations such as environmental, demographic, and other technical factors, and can be more severe when we combine different data sources.They have been shown to have substantial effects on gene expression levels and failing to taking them into consideration may lead to wrong scientific conclusions (Leek and Storey, 2007).When the data are aggregated from multiple sources, it remains an open problem on what is the best normalization practice.
Even with the systematic biases removed, another challenge is to conduct large-scale tests to pick important genes, proteins or SNPS.In testing the significance of thousands of genes, classical methods of controlling the probability of making one falsely discovered gene are no longer suitable and alternative procedures have been designed to control the false discovery rates (Benjamini and Hochberg, 1995;Storey, 2003;Schwartzman et al., 2009;Efron, 2010;Fan et al., 2012b) and to improve the power of the tests (Fan et al., 2012b).These technologies, though high-throughput in measuring the expression levels of tens of thousands of genes, remain low-throughput in surveying biological contexts (e.g., novel cell types, tissues, diseases, etc.).
An additional challenge in genomic data analysis is to model and explore the underlying heterogeneity of the aggregated datasets.Due to technology limitations and resource constraints, a single lab usually can only afford performing experiments for no more than a few cell types.This creates a major barrier for comprehensively characterizing gene regulation in all biological contexts, which is a fundamental goal of functional genomics.On the other hand, the NCBI Gene Expression Omnibus (GEO) and other public databases have cumulated more than 500, 000 gene expression profiles, including microarray, exon array, and RNA-seq samples from thousands of biological contexts.Public ChIP-chip and ChIP-seq data generated by different labs for different proteins and in different contexts are also steadily growing.Together, these public data contain enormous amounts of information that have not been fully exploited so far.Massive data aggregated from these public databases shed light on systematically studying many biological contexts in a highthroughput way.However, how to systematically explore the underlying heterogeneity and unveil the commonality across different sub-populations remains an active research area.

Neuroscience
Many important diseases, including Alzheimer's disease, Schizophrenia, Attention Deficit Hyperactive Disorder, Depression, and Anxiety, have been shown to be related to brain connectivity networks.Understanding the hierarchical, complex, functional network organization of the brain is a necessary first step to explore how the brain changes with disease.Rapid advances in neuroimaging techniques, such as functional magnetic resonance image (fMRI) and positron emission tomography (PET) as well as electrophysiology, provide great potential for the study of functional brain networks, i.e., the coherence of the activities among different brain regions (Jonides et al., 2006).
Take fMRI for example.It is a non-invasive technique for determining the neural correlates of mental processes in humans.During the past decade, this technique has become a leading method in the fields of cognitive and physiological neuroscience and kept producing massive amounts of high-resolution brain images.These images enable us to explore the association between brain connectivity and potential responses such as disease or psychological status.The fMRI data are massive and very high dimensional.Due to its non-invasive feature, everyday many fMRI machines keep scanning different subjects and constantly produce new imaging data.For each data point, the subject's brain is scanned for hundreds of times.Therefore, it is a 3D time-course image which contains more than hundreds of thousands of voxels.At the same time, the fMRI images are noisy due to its technological limit and possible head motion of the subjects.Analyzing such high dimensional and noisy data poses great challenges to statisticians and neuroscientists.
Similar to the field of genomics, an important Big Data problem in neuroscience is to aggregate datasets from multiple sources.Brain imaging data sharing is becoming more and more frequent nowadays (Visscher and Weissman, 2011).Primary sources of fMRI data arise from the International Data Sharing Initiative (INDI) and the 1,000 Functional Connectomes Project, Autism Brain Imaging Data Exchange (ABIDE) and ADHD-200 datasets.These international efforts have compiled thousands of resting-state fMRI scans along with complimentary structural scans.The largest of the datasets is the 1,000 Functional Connectomes Project, which focuses on healthy adults and includes limited covariate information on age, gender, handedness and image quality.The ADHD-200 dataset is similarly structured, yet includes diagnostic information on disease status such as human IQ.The ABIDE dataset is similar to the ADHD-200 dataset, with diagnostic autism and symptom severity information.However, it has a greater balance between diseased and non-diseased subjects.These large datasets pose great opportunities as well as new challenges.
One of the main challenges, as in the area of genomics, is to remove the systematic biases caused by experimental variations and data aggregations.Moreover, statisticallycontrolled inclusion of a subject in a group study, i.e., testing whether a person should be rejected as outlier data, is often poorly conducted (Fritsch et al., 2012) and voxels cannot be perfectly aligned across different experiments in different laboratories.Therefore, the collected data contain many outliers and missing values.These issues make data preprocessing and analysis significantly more complicated.Many traditional statistical procedures are not well-suited in this noisy high dimensional settings and new statistical thinking is crucially needed.

Economics and Finance
Over the past decade, more and more corporations are adopting the data-driven approach to conduct more targeted services, reduce risks, and improve performance.They are implementing specialized data analytics programs to collect, store, manage, and analyze large datasets from a range of sources to identify key business insights that can be exploited to support better decision making.For example, available financial data sources include stock prices, currency and derivative trades, transaction records, high-frequency trades, unstructured news and texts, consumers' confidence and business sentiments buried in social media and internet, among others.Analyzing these massive datasets helps measuring firms risks as well as systematic risks.It requires professionals who are familiar with sophisticated statistical techniques in portfolio management, securities regulation, proprietary trading, financial consulting, and risk management.
Analyzing large panel of economic and financial data is challenging.For example, as an important tool in analyzing the joint evolution of macroeconomics time series, the conventional vector autoregressive (VAR) model usually includes no more than ten variables, given the fact that the number of parameters grows quadratically with the size of the model.However, nowadays econometricians need to analyze multivariate time series with more than hundreds of variables.Incorporating all information into the VAR model will cause severe overfitting and bad prediction performance.One solution is to resort to sparsity assumptions, under which new statistical tools have been developed (Song and Bickel, 2011;Han and Liu, 2013b).
Another example is portfolio optimization and risk management (Cochrane, 2001;Dempster, 2002).In this problem, estimating the covariance and inverse covariance matrices of the returns of the assets in the portfolio plays an important role.Suppose that we have 1,000 stocks to be managed.There are 500, 500 covariance parameters to be estimated.Even if we could estimate each individual parameter accurately, the cumulated error of the whole matrix estimation can be large under matrix norms.This requires new statistical procedures.See, for example, Stock and Watson (2002); Bai and Ng (2002); Bai (2003); Forni et al. (2005); Fan et al. (2008); Bickel and Levina (2008); Cai and Liu (2011a); Agarwal et al. (2012b); Liu et al. (2012a); Xue and Zou (2012); Liu et al. (2012b); Fan et al. (2013); Pourahmadi (2013) on estimation large covariance matrices and their inverse.

Other Applications
Big Data have numerous other applications.Taking social network data analysis for example, massive amount of social network data are being produced by Twitter, Facebook, LinkedIn, and YouTube.These data reveal numerous individual's characteristics and have been exploited in various fields.For example, Aramaki et al. (2011) used these data to predict influenza epidemic; Bollen et al. (2011) used these data to predict the stock market trend; and Asur and Huberman (2010) used the social network data to predict box-office revenues for movies.In addition, the social media and internet contain massive amount of information on the consumer preferences and confidences, leading economics indicators, business cycles, political attitudes, and the economic and social states of a society.It is anticipated that the social network data will continue to explode and be exploited for many new applications.
Several other new applications that are becoming possible in the Big Data era include: (i) Personalized services: With more personal data collected, commercial enterprises are able to provide personalized services adapt to individual preferences.For example, Target (a retailing company in the United States) are able to predict a customer's need by analyzing the collected transaction records.
(ii) Internet security: When a network-based attack takes place, historical data on network traffic may allow us to efficiently identify the source and targets of the attack.
(iii) Personalized medicine: More and more health-related metrics such as individual's molecular characteristics, human activities, human habits, and environmental factors are now available.Using these pieces of information, it is possible to diagnose an individual's disease and select individualized treatments.
(iv) Digital humanities: Nowadays many archives are being digitized.For example, Google has scanned millions of books and identified about every word in every one of those books.This produces massive amount of data and enables addressing topics in the humanities, such as mapping the transportation system in ancient Roman, visualizing the economic connections of ancient China, studying how natural languages evolve over time, or analyzing historical events.

Salient Features of Big Data
Big Data create unique features that are not shared by the traditional datasets.These features pose significant challenges to data analysis and motivate the development of new statistical methods.Unlike traditional datasets where the sample size is typically larger than the dimension, Big Data are characterized by massive sample size and high dimensionality.First, we will discuss the impact of large sample size on understanding heterogeneity: On one hand, massive sample size allows us to unveil hidden patterns associated with small sub-populations and weak commonality across the whole population.On the other hand, modeling the intrinsic heterogeneity of Big Data requires more sophisticated statistical methods.Secondly, we discuss several unique phenomena associated with high dimensionality, including noise accumulation, spurious correlation, and incidental endogeneity.These unique features make traditional statistical procedures inappropriate.Unfortunately, most high-dimensional statistical techniques address only noise accumulation and spurious correlations issues, but not incidental endogeneity.They are based on exogeneity assumptions that often can not be validated by collected data, due to incidental endogeneity.

Heterogeneity
Big data are often created via aggregating many data sources corresponding to different sub-populations.Each sub-population might exhibit some unique features not shared by others.In classical settings where the sample size is small or moderate, data points from small sub-populations are generally categorized as "outliers" and it is hard to systematically model them due to insufficient observations.However, in the Big data era, the large sample size enables us to better understand heterogeneity, shedding lights toward studies such as exploring the association between certain covariates (e.g., genes or SNPs) and rare outcomes (e.g., rare diseases or diseases in small populations) and understanding why certain treatments (e.g.chemotherapy) benefit a subpopulation and harm another subpopulation.
To better illustrate this point, we introduce the following mixture model for the population: where λ j ≥ 0 represents the proportion of the j-th subpopulation, p j (y; θ j (x)) is the probability distribution of the response of the j-th subpopulation given the covariates x with θ j (x) as the parameter vector.In practice, many subpopulation are rarely observed, i.e., λ j is very small.When the sample size n is moderate, nλ j can be small, making it infeasible to infer the covariate-dependent parameters θ j (x) due to the lack of information.However, because Big Data are characterized by large sample size n, the sample size nλ j for the j-th subpopulation can be moderately large even if λ j is very small.This enables us to more accurately infer about the sub-population parameters θ j (•).In short, the main advantage brought by Big Data is to understand heterogeneity of sub-populations such as the benefits of certain personalized treatments, which are infeasible when sample size is small or moderate.Big Data also allow us to unveil weak commonality across whole population, thanks to large sample sizes.For example, the benefit of one drink of red wine per night on heart can be difficult to assess without large sample.Similarly, health risks to exposure of certain environmental factors can only be more convincingly evaluated when the sample sizes are sufficiently large.
Besides the aforementioned advantages, the heterogeneity of Big Data also poses significant challenges to statistical inference.Inferring the mixture model in (3.1) for large datasets requires sophisticated statistical and computational methods.In low dimensions, standard techniques such as the expectation-maximization algorithm for finite mixture models can be applied.In high dimensions, however, we need to carefully regularize the estimating procedure to avoid overfitting or noise accumulation and to devise good computation algorithms (Khalili and Chen, 2007;Städler et al., 2010)

Noise Accumulation
Analyzing Big Data requires us to simultaneously estimate or test many parameters.These estimate errors accumulate when a decision or prediction rule depends on a large number of such parameters.Such a noise accumulation effect is especially severe in high dimensions and may even dominate the true signals.It is usually handled by the sparsity assumption (Donoho, 2000;Hastie et al., 2009;Bühlmann and van der Geer, 2011).
Take high-dimensional classification for instance.Poor classification is due to the existence of many weak features that do not contribute to the reduction of classification error (Fan and Fan, 2008).As an example, we consider a classification problem where the data come from two classes: (3.2) We want to construct a classification rule which classifies a new observation Z ∈ R d into either the first or the second class.To illustrate the impact of noise accumulation in classification, we set n = 100 and d = 1, 000.We set µ 1 = 0 and µ 2 to be sparse, i.e., only the first 10 entries of µ 2 is nonzero with value 3, all the other entries are zero.
Figure 1 plots the first two principal components by using the first m = 2, 40, 200 features and the whole 1,000 features.As illustrated in these plots, when m = 2 we get high discriminative power.However, the discriminative power becomes very low when m is too large due to noise accumulation.The first 10 features contribute to classifications and the remaining features do not.Therefore, when m > 10, procedures do not get any additional signals, but accumulate noises: The larger m, the more noise accumulation, which makes the classification procedure deteriorates with dimensionality.For m = 40, the accumulated signals compensate the accumulated noise so that the first two-principal components still have good discriminative power.When m = 200, the accumulated noise exceeds the signal gains.
The above discussion motivates the usage of sparse models and variable selection to overcome the effect of noise accumulation.For example, in the classification model (3.2), instead of using all the features, we could select a subset of features which attain the best signal-to-noise ratio.Such a sparse model provides more improved classification performance (Hastie et al., 2009;Bühlmann and van der Geer, 2011).In other words, variable (a) m = 2 (b) m = 40 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −2 0 2 4 6 −2 0 2 4 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −4 −2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −5.0 −2.5 0.0 2.5 5.0 −6 −4 −2 0 2 4 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −10 −5 selection plays a pivotal role in overcoming noise accumulation in classification and regression prediction.However, variable selection in high dimensions is challenging due to spurious correlation, incidental endorgeneity, heterogeneity, and measurement errors.

Spurious Correlation
High dimensionality also brings spurious correlation, referring to the fact that many uncorrelated random variables may have high sample correlations in high dimensions.Spurious correlation may casue false scientific discoveries and wrong statistical inferences.
Consider the problem of estimating the coefficient vector β of a linear model where y ∈ R n represents the response vector, X = [x 1 , . . ., x n ] T ∈ R n×d represents the design matrix, ∈ R n represents an independent random noise vector, and I d is the d by d identity matrix.To cope with the noise accumulation issue, when the dimension d is comparable to or larger than the sample size n, it is popular to assume that only a small number of variables contribute to the response, i.e., β is a sparse vector.Under this sparsity assumption, variable selection can be conducted to avoid noise accumulation, improve the performance of prediction, as well as enhance the interpretability of the model with parsimonious representation.
In high dimensions, even for a model as simple as (3.3), variable selection is challenging due to the presence of spurious correlation.In particular, Fan and Lv (2008) showed that, when the dimensionality is high, the important variables can be highly correlated to several spurious variables which are scientifically unrelated.We consider a simple example to illustrate this phenomenon.Let x 1 , . . ., x n be n independent observations of a d-dimensional Gaussian random vector X = (X 1 , . . ., X d ) T ∼ N d (0, I d ).We repeatedly simulate the data with n = 60 and d = 800 and 6, 400 for 1,000 times.Figure 2 (a) shows the empirical distribution of the maximum absolute sample correlation coefficient between the first variable with the remaining ones defined as where Corr (X 1 , X j ) is the sample correlation between the variables X 1 and X j .We see that the maximum absolute sample correlation becomes higher as dimensionality increases.Furthermore, we can compute the maximum absolute multiple correlation between X 1 and linear combinations of several irrelevant spurious variables: (3.5) Using the same configuration as in Figure 2 (a), Figure 2 (b) plots the empirical distribution of the maximum absolute sample correlation coefficient between X 1 and j∈S β j X j , where S is any size 4 subset of {2, . . ., d} and β j is the least squares regression coefficient of X j when regressing X 1 on {X j } j∈S .Again, we see that even though X 1 is utterly independent of X 2 , . . ., X d , the correlation between X 1 and the closest linear combination of any four variables of {X j } j =1 to X 1 can be very high.We refer to Cai and Jiang (2012) and Fan et al. (2012a) about more theoretical results on characterizing the orders of r.
The spurious correlation has significant impact on variable selection and may lead to false scientific discoveries.Let X S = (X j ) j∈S be the sub-random vector indexed by S and S the selected set that has the higher spurious correlation with X 1 as in Figure 2.For example, when n = 60 and d = 6, 400, we see that X 1 is practically indistinguishable from X S for a set S with | S| = 4.If X 1 represents the expression level of a gene that is Here the dimension d is 800 and 6,400, the sample size n is 60.The result is based on 1,000 simulations.
responsible for a disease, we can not distinct it from the other 4 genes in S that have a similar predictive power although they are scientifically irrelevant.
Besides variable selection, spurious correlation may also lead to wrong statistical inference.We explain this by considering again the same linear model as in (3.3).Here we would like to estimate the standard error σ of the residual, which is prominently featured in statistical inferences of regression coefficients, model selection, goodness-of-fit test, and marginal regression.Let S be a set of selected variables and P S be the projection matrix on the column space of X S .The standard residual variance estimator, based on the selected variables, is The estimator (3.6) is unbiased when the variables are not selected by data and the model is correct.However, the situation is completely different when the variables are selected based on data.In particular, Fan et al. (2012a) showed that when there are many spurious variables, σ 2 is seriously underestimated, which leads further to wrong statistical inferences including model selection or significance tests, and false scientific discoveries such as finding wrong genes for molecular mechanisms.They also propose a refitted cross-validation method to attenuate the problem.

Incidental Endogeneity
Incidental endogeneity is another subtle issue raised by high dimensionality.In a regression setting Y = d j=1 β j X j + ε, the term "endogeneity" (Engle et al., 1983) means that some predictors {X j } correlate with the residual noise ε.The conventional sparse model assumes Y = j β j X j + ε, and E(εX j ) = 0 for j = 1, . . ., d, (3.7) with a small set S = {j : β j = 0}.The exogenous assumption in (3.7) that the residual noise ε is uncorrelated with all the predictors is crucial for validity of most existing statistical procedures, including variable selection consistency.Though this assumption looks innocent, it is easy to be violated in high dimensions as some of variables {X j } are incidentally correlated with ε, making most high-dimensional procedures statistically invalid.
To explain the endogeneity problem in more details, suppose that unknown to us, the response Y is related to three covariates as follows: with EεX j = 0, for j = 1, 2, 3.
In the data collection stage, we do not know the true model, and therefore collect as many covariates that are potentially related to Y as possible, in hope to include all members in S in (3.7).Incidentally, some of those X j 's (for j = 1, 2, 3) might be correlated with the residual noise ε.This invalidates the exogenous modeling assumption in (3.7).In fact, the more covariates are collected or measured, the harder this assumption is satisfied.Unlike spurious correlation, incidental endogeneity refers to the genuine existence of correlations between variables unintentionally, both due to high-dimensionality.The former is analogous to find two persons look alike but have no genetic relation, whereas the latter is similar to bumping into an acquaintance, both easily occurring in a big city.More generally, endogeneity happens as a result of selection biases, measurement errors, and omitted variables.These phenomena arise frequently in the analysis of Big Data, mainly due to two reasons: • With the benefit of new high-throughput measurement techniques, scientists are able to and tend to collect as many features as possible.This accordingly increases the possibility that some of them might be correlated to the residual noise, incidentally.
• Big Data are usually aggregated from multiple sources with potentially different data generating schemes.This increases the possibility of selection bias and measurement errors, which also cause potential incidental endogeneity.
Whether incidental endogeneity appears in real datasets and how shall we test it in practice?We consider a genomics study in which 148 microarray samples are downloaded from GEO database and ArrayExpress.These samples are created under Affymetrix HGU133a Here ε represents the residual noise after the Lasso fit.We provide the distributions of the sample correlations using both the raw data and permuted data.
platform for human subjects with prostate cancer.The obtained dataset contains 22,283 probes, corresponding to 12,719 genes.In this example, we are interested in the gene named "Discoidin domain receptor family, member 1" (abbreviated as DDR1).DDR1 encodes receptor tyrosine kinases, which plays an important role in communication of cells with their microenvironment.DDR1 is known to be highly related to the prostate cancer (Valiathan et al., 2012) and we wish to study its association with other genes in patients with prostate cancer.We took the gene expressions of DDR1 as the response variable Y and the expressions of all the remaining 12,718 genes as predictors.The left panel of Figure 3 draws the empirical distribution of the correlations between the response and individual predictors.
To illustrate the existence of endogeneity, we fit an L 1 -penalized least squares regression (Lasso) on the data and the penalty is automatically selected via ten-fold cross validation (37 genes are selected).We then refit an ordinary least squares regression on the selected model to calculate the residual vector.In the right panel of Figure 3, we plot the empirical distribution of the correlations between the predictors and the residuals.We see the residual noise is highly correlated with many predictors.To make sure these correlations are not purely caused by spurious correlation, we introduce a "null distribution" of the spurious correlations by randomly permuting the orders of rows in the design matrix such that the predictors are indeed independent of the residual noise.By comparing the two distributions, we see that the distribution of correlations between predictors and residual noise on the raw data (labeled "raw data") has heaviers tail than that on the permuted data (labeled "permuted data").This result provides stark evidence of endogeneity in this data.
The above discussion shows that incidental endogeneity is likely to be present in Big Data.The problem of dealing with endogenous variables is not well understood in high dimensional statistics.What is the consequence of this endogeneity?Fan and Liao (2012) showed that endogeneity causes model selection inconsistency.In particular, they provided thorough analysis to illustrate the impact of endogeneity on high dimensional statistical inference and proposed alternative methods to conduct linear regression with consistency guarantees under weaker conditions.See also the next section.

Impact on Statistical Thinking
As has been shown in the previous section, massive sample size and high dimensionality bring heterogeneity, noise accumulation, spurious correlation, and incidental endogeneity.These features of Big Data make traditional statistical methods invalid.In this section we introduce new statistical methods that can handle these challenges.For an overview, see Hastie et al. (2009) and Bühlmann and van der Geer (2011).

Penalized Quasi-likelihood
To handle noise accumulation issue, we assume the model parameter β as in (3.3) is sparse.The classical model selection theory, according to Akaike (1974), suggests to choose a parameter vector β that minimizes negative penalized quasi-likelihood: where QL(β) is the quasi-likelihood of β and • 0 represents the L 0 pseudo-norm (i.e., the number of nonzero entries in a vector).Here λ > 0 is a regularization parameter that controls the bias-variance tradeoff.The solution to the optimization problem in (4.1) has nice statistical properties (Barron et al., 1999).However, it is essentially combinatoric optimization and does not scale to large-scale problems.The estimator in (4.1) can be extended to a more general form: where the term n (β) measures the goodness of fit of the model with parameter β and d j=1 P λ,γ (β j ) is a sparsity-inducing penalty that encourages sparsity, in which λ is again the tuning parameter that controls the bias-variance tradeoff and γ is a possible finetune parameter which controls the degree of concavity of the penalty function (Fan and Li, 2001).Popular choices of the penalty function P λ,γ (•) include the hard-thresholding penalty (Antoniadis, 1997;Antoniadis and Fan, 2001), soft-thresholding penalty (Donoho  and Johnstone, 1994;Tibshirani, 1996), smoothly clipped absolution deviation (SCAD, Fan and Li ( 2001)), and minimax concavity penalty (MCP, Zhang (2010)).Figure 4 visualizes these penalty functions for λ = 1.We see that all penalty functions are folded concave, but the soft-thresholding (L 1 -) penalty is also convex.The parameter γ in SCAD and MCP controls the degree of concavity.From Figure 4(c) and Figure 4(d), we see that a smaller value of γ results in more concave penalties.When γ becomes larger, SCAD and MCP converge to the soft-thresholding penalty.MCP is a generalization of the hard-thresholding penalty which corresponds to γ = 1.How shall we choose among these penalty functions?In applications, we recommend to use either SCAD or MCP thresholding since they combine the advantages of both hardand soft-thresholding operators.Many efficient algorithms have been proposed for solving the optimization problem in (4.2) with the above four penalties.See Section 5.

Sparsest Solution in High Confidence Set
The penalized quasi-likelihood estimator (4.2) is somewhat mysterious.A closely related method is the the sparsest solution in high confidence set, introduced in the recent book chapter by Fan (2013), which has much better statistical intuition.It is a generally applicable principle that separates the data information and the sparsity assumption.
Suppose that the data information is summarized by the function n (β) in (4.2).This can be a likelihood, quasi-likelihood, or loss function.The underlying parameter vector β 0 usually satisfies (β 0 ) = 0, where (•) is the gradient vector of the expected loss function (β) = E n (β).Thus, a natural confidence set for β 0 is where • ∞ is the L ∞ -norm of a vector and γ n is chosen so that we have confidence level at least 1 − δ n , namely The confidence set C n is called high-confidence set since δ n → 0. In theory, we can take any norm in constructing the high-confidence set.We opt for the L ∞ norm, as it produces a convex confidence set C n when n (•) is convex.The high-confidence set is a summary of the information we have for the parameter vector β 0 .It is not informative in high-dimensional space.Take, for example, the linear model (3.3) with the quadratic loss n (β) = y − Xβ 2 2 .The high-confidence set is then where we take γ n ≥ X T ε ∞ so that δ n = 0.If in addition β 0 is assumed to be sparse, then a natural solution is the intersection of these two pieces of information, namely, finding the sparsest solution in the high-confidence set: This is a convex optimization problem when (•) is convex.For the linear model with the quadratic loss, it reduces to the Dantzig selector (Candes and Tao, 2007).
There are many flexibilities in defining the sparsest solution in high-confidence set.First of all, we have a choice of the loss function n (•).We can regard n (β) = 0 as the estimation equations (Liang and Zeger, 1986) and define directly the high confidence set (4.3) from the estimation equations.Secondly, we have many ways to measure the sparsity.For example, we can use a weighted L 1 -norm to measure the sparsity of β in (4.5).By proper choices of estimating equations in (4.3) and measure of sparsity in (4.5), Fan (2013) showed that many useful procedures can be regarded as the sparsest solution in the high-confidence set.
For example, CLIME (Cai et al., 2011) for estimating sparse precision matrix in Gaussian graphic model and the linear programming discriminant rule (Cai and Liu, 2011b) for sparse high-dimensional classification are both the sparsest solution in high confidence set.Fan (2013) also provided a general convergence theory for such a procedure under a condition similar to the restricted eigenvalue condition in Bickel et al. (2009b).Finally, the idea is applicable to the problems with measurement errors or even endogeneity.In this case, the high-confidence set will be defined accordingly to accommodate the measurement errors or endogeneity.See, for example, Gautier and Tsybakov (2011).

Independence Screening
An effective variable screening technique based on marginal screening has been proposed by Fan and Lv (2008).They aim at handling ultra-high dimensional data for which the aforementioned penalized quasi-likelihood estimators become computationally infeasible.For such cases, Fan and Lv (2008) proposed to first use marginal regression to screen variables, reducing the original large-scale problem to a moderate-scale statistical problem so that more sophisticated methods for variable selection can be applied.The proposed method, named sure independence screening, is computationally very attractive.It has been shown to possess sure screening property and to have some theoretical advantages over Lasso (Fan and Song, 2010;Genovese et al., 2012) There are two main ideas of sure independent screening: (i) It uses the marginal contribution of a covariate to probe its importance in the joint model; (ii) Instead of selecting the most important variables, it aims at removing variables that are not important.For example, assuming each covariate has been standardized, we denote β M j the estimated regression coefficient in a univariate regression model.The set of covariates that survive the marginal screening is defined as for a given threshold δ.One can also measure the importance of a covariate X j by using its deviance reduction.For the least-squares problem, both methods reduce to ranking importance of predictors by using the magnitudes of their marginal correlations with the response Y .Fan and Lv (2008) and Fan and Song (2010) gave conditions under which sure screening property can be established and false selection rates are controlled.
Since the computational complexity of sure screening scales linearly with the problem size, the idea of sure screening is very effective in dramatic reduction of the computational burden of Big Data analysis.It has been extended in various directions.For example, generalized correlation screening was used in Hall and Miller (2009), nonparametric screening was proposed by Fan et al. (2011), and principled sure independence screening was introduced in Zhao and Li (2012).In addition, Li et al. (2012b) utilized the distance correlation to conduct screening, Li et al. (2012a) employed rank correlation, and Fan et al. (2009) proposed an iteratively screening and selection method.
Independent screening has never examined multivariate effect of variables on the response variable nor has it used the covariance matrix of variables.An extension of this is to use multivariate screening, which examines the contributions of small groups of variables together.This allows us to examine the synergy of small groups of variables to the response variable.However, the bivariate screening already involves O(d 2 ) submodels, which can be prohibitive in computation.Covariance assist screening and estimation in Ke et al. (2012) can be adapted here to prevent examining all bivariate or multivariate submodels.Another possible extension is to develop conditional screening techniques, which rank variables according to their conditional contributions given a set of variables.

Dealing with Incidental Endogeneity
Big Data is prone to incidental endogeneity that makes most popular regularization method invalid.It is accordingly important to develop methods that can handle endogeneity in high dimensions.More specifically, let's consider the high dimensional linear regression model (3.7).Fan and Liao (2012) showed that for any penalized estimators to be variable selection consistent, a necessary condition is E(εX j ) = 0 for j = 1, . . ., d. (4.7) As discussed in the Section 3, the condition in (4.7) is too restrictive for real-world applications.Letting S = {j : β j = 0} be the set of important variables, with non-vanishing components in β, a more realistic model assumption should be In the paper by Fan and Liao (2012), they considered an even weaker version of Equation (4.8), called the "over identification" condition, such as EεX j = 0 and EεX 2 j = 0 for j ∈ S. (4.9) Under condition (4.9), Fan and Liao (2012) showed that the classical penalized least squares methods such as lasso, SCAD, and MCP, are no longer consistent.Instead, they introduced the Focused Generalized Methods of Moments (FGMM) by utilizing the over identification conditions and proved that the FGMM consistently selects the set of variables S. We do not go into the technical details here but illustrate this by an example.We continue to explore the gene expression data in Section 3.4.We again treat gene DDR1 as response and other genes as predictors, and apply the FGMM instead of Lasso.By cross validation, the FGMM selects 18 genes.The left panel of Figure 5 shows the distribution of the sample correlations between the genes X j (j = 1, . . ., 12718) and the residuals ε after the FGMM fit.Here we find that many correlations are non-zero, but it does not matter, because we require only (4.9).To verify (4.9), the right panel of Figure 5 shows the distribution of the sample correlations between the 18 selected genes (and their squares) and the residuals.The sample correlations between the selected genes and residuals are zero, and the sample correlations between the squared covariates and residuals are small.Therefore, the modeling assumption is consistent to our model diagnostics.

Impact on Computing Infrastructure
The massive sample size of Big Data fundamentally challenges the traditional computing infrastructure.In many applications we need to analyze internet-scale data which contain billions or even trillions of data points, which even makes a linear pass of the whole dataset unaffordable.In addition, such data could be highly dynamic and infeasible to be stored in a centralized database.The fundamental approach to store and process such data is to divide and conquer.The idea is to partition a large problem into more tractable and independent subproblems.Each subproblem is tackled in parallel by different processing units.Intermediate results from each individual worker are then combined to yield the final output.In small scale, such divide-and-conquer paradigm can be implemented either by multi-core computing or grid computing.However, in very large scale, it poses fundamen-tal challenges to computing infrastructure.For example, when millions of computers are connected to scale out to large computing tasks, it is quite likely some computers may die during the computing.In addition, given a large computing task, we want to distribute it evenly to many computers and make the workload balanced.Designing very large scale, high adaptive and fault-tolerant computing systems is extremely challenging and motivates the outcome of new and reliable computing infrastructure that supports massively parallel data storage and processing.In this section we take Hadoop as an example to introduce basic software and programming infrastructure for Big Data processing.• HDFS is a self-healing, high-bandwidth, clustered storage file system, and • MapReduce is a distributed programming model developed by Google.
We dart with explaining HDFS and MapReduce in the next two subsections.Besides these two key components, a typical Hadoop release contains many other components.For example, as is shown in Figure 6, the Cloudera's open-source Hadoop distribution also includes HBase, Hive, Pig, Oozie, Flume and Sqoop.More details about these extra components are provided in the online Cloudera technical documents.After introducing the Hadoop, we also briefly explain the concepts of cloud computing in Sections 5.3.

Hadoop Distributed File System
HDFS is a distributed file system designed to host and provide high-throughput access to large datasets which are redundantly stored across multiple machines.In particular, it ensures Big Data's durability to failure and high availability to parallel applications.
As a motivating application, suppose we have a large data file containing billions of records and we want to query this file frequently.If many queries are submitted simultaneously (e.g., the Google search engine), the usual file system is not suitable due to the I/O limit.HDFS solves this problem by dividing a large file into small blocks and store them in different machines.Each machine is called a DataNode.Unlike most block-structured file systems which use a block size on the order of 4 or 8 KB, the default block size in HDFS is 64MB, which allows HDFS to reduce the amount of metadata storage required per file.Furthermore, HDFS allows for fast streaming reads of data by keeping large amounts of data sequentially laid out on the hard disk.The main tradeoff of this decision is that HDFS expects the data to be read sequentially (instead of being read in a random access fashion).
The data in HDFS can be accessed via a "write once and read many" approach: The metadata structures (e.g., the file names and directories) are allowed to be simultaneously modified by many clients.It is important that this meta information is always synchronized and stored reliably.All the metadata is maintained by a single machine, called the NameNode.Because of the relatively low amount of metadata per file (it only tracks file names, permissions, and the locations of each block of each file), all such information can be stored in the main memory of the NameNode machine, allowing fast access to the metadata.An illustration of the whole HDFS architecture is provided in Figure 7.To access or manipulate a data file, a client contacts the NameNode and retrieves a list of locations for the blocks that comprise the file.These locations identify the DataNodes which hold each block.Clients then read file data directly from the DataNode servers, possibly in parallel.The NameNode is not directly involved in this bulk data transfer, keeping its working load to a minimum.HDFS has a built-in redundancy and replication feature which secures that any failure of individual machines can be recovered without any loss of data (e.g., each DataNode has 3 copies by default).The HDFS automatically balances its load whenever a new DataNode is added to the cluster.We also need to safely store the NameNode information by creating multiple redundant systems, which allows the important metadata of the file system be recovered even if the NameNode itself crashes.

MapReduce
MapReduce is a programming model for processing large datasets in a parallel fashion.We use an example to explain how MapReduce works.Suppose we are given a symbol sequence (e.g., "ATGCCAATCGATGGGACTCC"), and the task is to write a program that counts the number of each symbol.The simplest idea is to read a symbol, add it into a hash table with key as the symbol and set value to its number of occurrences.If the symbol is not in the hash table yet, then add the symbol as a new key to the hash and set the corresponding value to 1.If the symbol is already in the hash table, then increase the value by 1.This program runs in a serial fashion and the time complexity scales linearly with the length of the symbol sequence.Everything looks simple so far.However, imagine if instead of a simple sequence, we need to count the number of symbols in the whole genomes of many biological subjects.Serial processing of such a huge amount of information is time consuming.So, the question is how can we use parallel processing units speed up the computation.
The idea of MapReduce is illustrated in Figure 8.We initially split the original sequence into several files (e.g., 2 files in this case).We further split each file into several subsequences (e.g., 2 subsequences in this case) and "map" the number of each symbol in each subsequence.The outputs of the mapper are (key, value)-pairs.We then gather together all output pairs of the mappers with the same key.Finally, we use a "reduce" function to combine the values for each key.This gives the desired output: The Hadoop MapReduce contains three stages, which are listed as follows.

First Stage: Mappling
The first stage of a MapReduce program is called mapping.In this stage, a list of data elements are provided to a "mapper" function to be transformed into (key, value)-pairs.For example, in the above symbol counting problem, the mapper function simply transforms

Intermediate Stages: Shuffling and Sorting
After the mapping stage, the program exchanges the intermediate outputs from the mapping stage to different "reducers".This process is called shuffling.A different subset of the intermediate key space is assigned to each reduce node.These subsets (known as "partitions") are the inputs to the next reducing step.Each map task may send (key, value)-pairs to any partition.All pairs with the same key are always grouped together on the same reducer regardless of which mappers they are coming from.Each reducer may process several sets of pairs with different keys.In this case, different keys on a single node are automatically sorted before they are fed into the next Reducing step.

Final Stage: Reducing
In the final reducing stage, an instance of user-provided code is called for each key in the partition assigned to a reducer.The inputs are a key and an iterator over all the values associated with the key.These values returned by the iterator could be in an undefined order.In particular, we have one output file per executed reduce task.The Hadoop MapReduce builds on the HDFS and inherits all the fault-tolerance properties of HDFS.In general, Hadoop is deployed on very large-scale clusters.One example is shown in Figure 9.

Cloud Computing
Cloud computing revolutionizes modern computing paradigm.It allows everything -from hardware resources, software infrastructure to datasets -to be delivered to data analysts as a service wherever and whenever needed.Figure 10 illustrates different building components of cloud computing.The most striking feature of cloud computing is its elasticity and ability to scale up and down, which makes it suitable for storing and processing Big Data.

Impact on Computational Methods
Big data are massive and very high dimensional, which pose significant challenges on computing and paradigm shifts on large-scale optimization (Boyd and Vandenberghe, 2004;Boyd et al., 2011).On one hand, directly application of penalized quasi-likelihood estimators on high dimensional data requires us to solve very large-scale optimization problems.Optimization with a large amount of variables is not only expensive but also suffers from slow numerical rates of convergence and instability.Such a large-scale optimization is generally regarded as a mean, not the goal of Big Data analysis.Scalable implementations of calculating the gradient of the objective function at each point.However, calculating gradient can be very time consuming when the dimensionality is high.Instead, Friedman et al. (2007) proposed to calculate the penalized pseudo-likelihood estimator using the pathwise coordinate descent algorithm, which can be viewed as a special case of the gradient descent algorithm: Instead of optimizing along the direction of the full gradient, it only calculates the gradient direction along one coordinate at each time.A beautiful feature of this is that even though the whole optimization problem does not have a closed-form solution, there exist simple closed-form solutions to all the univariate subproblems.The coordinate descent is computationally easy and has similar numerical convergence properties as gradient descent (Nesterov, 2012).Alternative first-order algorithms to coordinate descent have also been proposed and widely used, resulting in iterative shrinkage-thresholding algorithms (Daubechies et al., 2004;Beck and Teboulle, 2009).Prior to the coordinate descent algorithm, Efron et al. (2004) proposed the least angle regression (LARS) algorithm to the L 1 -penalized least squares problem.
When the penalty function P λ,γ (•) is nonconvex (e.g., SCAD and MCP), the objective function in (4.2) is no longer concave.Many algorithms have been proposed to solve this optimization problem.For example, Fan and Li (2001) proposed a local quadratic approximation (LQA) algorithm for optimizing nonconcave penalized likelihood.Their idea is to approximate the penalty term piece by piece using a quadratic function, which can be thought as a convex relaxation (majorization) to the nonconcave object function.With the quadratic approximation, a closed-form solution can be obtained.This idea is further improved by using a linear instead of a quadratic function to approximate the penalty term and leads to the local linear approximation (LLA) algorithm (Zou and Li, 2008).More specifically, given current estimate d ) T at the k th iteration for problem (4.2), by Taylor's expansion, (6.1) Thus, at the (k + 1) th iteration, we solve min where w k,j = P λ,γ (β (k) j ).Note that problem (6.2) is convex so that a convex solver can be used.Fan et al. (2008) suggested using initial values β (0) = 0, which corresponds to the unweighted L 1 penalty.This algorithm shares a very similar idea as in Candes et al. (2008), both of which can be regarded as implementations of minimization of the folded-concave penalized quasi-likelihood (Fan and Li, 2001) problem (4.2).If one further approximates the goodness of fit measure n (β) in (6.2) by a quadratic function via Taylor expansion, then the LARS algorithm (Efron et al., 2004) and pathwise coordinate descent algorithm (Friedman et al., 2007) can be used.
For the more general settings where the loss function n (•) may not be concave, Wang et al. (2013) proposed an approximate regularization path following algorithm for solving the optimization problem in (4.2).By integrating statistical analysis with computational algorithms, they provided explicit statistical and computational rates of convergence of any local solution obtained by the algorithm.Computationally, the approximate regularization path following algorithm attains a global geometric rate of convergence for calculating the full regularization path, which is fastest possible among all first-order algorithms in terms of iteration complexity.Statistically, they show that any local solution obtained by the algorithm attains the oracle properties with the optimal rates of convergence.The idea on studying statistical properties based on computational algorithms, which combine both computational and statistical analysis, represents an interesting future direction for Big Data.We also refer to Agarwal et al. (2012a) and Loh and Wainwright (2013) for researches in this direction.

Dimension Reduction and Random Projection
We introduce several dimension (data) reduction procedures in this section.Why do we need dimension reduction?Let's consider a dataset represented as a n by d real-value matrix D, which encodes information about n observations of d variables.In the Big data era, it is in general computationally intractable to directly make inference on the raw data matrix.Therefore, an important data preprocessing procedure is to conduct dimension reduction which finds a compressed representation of D that is of lower dimensions but preserves as much information in D as possible.
Principal component analysis is the most well known dimension reduction method.It aims at projecting the data onto a low dimensional orthogonal subspace that captures as much of the data variation as possible.Empirically, it calculates the leading eigenvectors of the sample covariance matrix to form a subspace U k ∈ R d×k .We then project the n×d data matrix D to this linear subspace to get a n×k data matrix D U k .This procedure is optimal among all the linear projection methods in minimizing the squared error introduced by the projection.However, conducting the eigenspace decomposition on the sample covariance matrix is computational challenging when both n and d are large: The computational complexity of PCA is O(d 2 n + d 3 ) (Golub and Van Loan, 2012), which is infeasible for very large datasets.
To handle the computational challenge raised by massive and high dimensional datasets, we need to develop methods that preserve the data structure as much as possible and is computational efficient for handling high dimensionality.Random projection (Johnson and Lindenstrauss, 1984) is an efficient dimension reduction technique for this purpose, and is closely related to the celebrated idea of compress sensing (Donoho, 2006;Tsaig and Donoho, 2006;Lustig et al., 2007;Figueiredo et al., 2007;Candès and Wakin, 2008).More specifically, random projection aims at finding a k-dimensional subspace of D such that the distances between all pairs of data points are approximately preserved.It achieves this goal by projecting the original data D onto a k-dimensional subspace using a random projection matrix with unit column norms.More specifically, let R ∈ R d×k be a random matrix with all the column Euclidean norms equal to one.We reduce the dimensionality of D from d to k by calculating matrix multiplication This procedure is very simple and the computational complexity of the random projection procedure is of order O(ndk), which scales only linearly with the problem size.
Theoretical justifications of random projection are based on two results.Johnson and Lindenstrauss (1984) showed that if points in a vector space are projected onto a randomly selected subspace of suitable dimensions, then the distance between the points are approximately preserved.This justifies the random projection when R is indeed a projection matrix.However, enforcing R to be orthogonal requires Gram-Schmidt algorithm, which is computationally expensive.In practice, Marks and Zurada (1994) showed that in high dimensions we do not need to enforce the matrix to be orthogonal.In fact, any finite number of high dimensional random vectors are almost orthogonal to each other.This result guarantees that R T R can be sufficiently close to the identity matrix.Achlioptas (2001) further simplified the random projection procedure by removing the unit column length constraint.
To illustrate the usefulness of random projection, we use the gene expression data in Section 3.4 to compare the performance of PCA and random projection in preserving the relative distances between pairwise data points.We extract the top 100, 500, and 2,500 genes with the highest marginal standard deviations, then apply PCA and random projection (RP) to reduce the dimensionality of the raw data to a small number k. Figure 11 shows the median errors in the distance between members across all pairs of data vectors.We see that, when dimensionality increases, random projections have more and more advantages over PCA in preserving the distances between sample pairs.One thing to note is that random projection is not the "optimal" procedure for traditional small scale problems.Accordingly, the popularity of this dimension reduction procedure indicates a new understanding of Big data.To balance the statistical accuracy and computational complexity, the suboptimal procedures in small or medium scale problems can be "optimal" in large scale.Moreover, the theory of random projection depends on the high dimensionality feature of Big data.This can be viewed as a blessing of dimensionality.
Besides PCA and random projection, there are many other dimension reduction methods, including latent semantic indexing (LSI) (Deerwester et al., 1990), discrete cosine transform (Rao et al., 2007), and CUR decomposition (Mahoney and Drineas, 2009).These methods have been widely used in analyzing large text and image datasets.Here "RP" stands for the random projection and "PCA" stands for the principal component analysis.

Conclusions and Future Perspectives
This paper discusses statistical and computational aspects of Big Data analysis.We selectively overview several unique features brought by Big Data and discuss some solutions.
Besides the challenge of massive sample size and high dimensionality, there are several other important features of Big Data worth equal attention.These include (1) Complex data challenge: Due to the fact that Big Data are in general aggregated from multiple sources, they sometime exhibit heavy tail behaviors with nontrivial tail dependence.
(2) Noisy data challenge: Big Data usually contain various types of measurement errors, outliers, and missing values.
(3) Dependent data challenge: In various types of modern data, such as financial time series, fMRI and time course microarray data, the samples are dependent with relatively weak signals.

Figure 1 :
Figure 1: Scatter plots of projections of the observed data (n = 100 from each class) onto the first two principal components of the best m-dimensional selected feature space.A projected data with "•" indicates the first class and " " indicates the second class.

Figure 2 :
Figure 2: Illustration of spurious correlation.(a): Distribution of the maximum absolute sample correlation coefficients between X 1 and {X j } j =1 .(b): Distribution of the maximum absolute sample correlation coefficients between X 1 and the closest linear projections of any 4 members of {Xj } j =1 to X 1 .Here the dimension d is 800 and 6,400, the sample size n is 60.The result is based on 1,000 simulations.

Figure 3 :
Figure 3: Illustration of incidental endogeneity on a microarry gene expression data.Left panel: Distribution of the sample correlation Corr(X j , Y ) (j = 1, • • • , 12, 718).Right panel: Distribution of the sample correlation Corr(X j , ε).Here ε represents the residual noise after the Lasso fit.We provide the distributions of the sample correlations using both the raw data and permuted data.

Figure 4 :
Figure 4: Visualization of the penalty functions.In all cases, λ = 1.For SCAD and MCP, different values of γ are chosen as shown in graphs.

Figure 5 :
Figure 5: Diagnostics of the modeling assumptions of the FGMM on a microarry gene expression data.Left panel: Distribution of the sample correlations Corr(X j , ε) (j = 1, • • • , 12, 718).Right panel: Distribution of the sample correlations Corr(X j , ε) and Corr(X 2 j , ε) for only 18 selected genes.Here ε represents the residual noise after the FGMM fit.

Figure 6 :
Figure 6: An illustration of Cloudera's open-source Hadoop distribution (source: cloudera website).Hadoop is a Java-based software framework for distributed data management and processing.It contains a set of open source libraries for distributed computing using the MapReduce programming model and its own distributed file system called HDFS.Hadoop automatically facilitates scalability and takes cares of detecting and handling failures.Core Hadoop has two key components: Core Hadoop = Hadoop distributed file system (HDFS) + MapReduce

Figure 7 :
Figure 7: An illustration of the HDFS architecture.

Figure 8 :
Figure 8: An illustration of the MapReduce paradigm for the symbol counting task.Mappers are applied to every element of the input sequences and emit intermediate (key, value)pairs.Reducers are applied to all values associated with the same key.Between the map and reduce stages are some intermediate steps involving distributed sorting and grouping.each symbol into the pair (symbol, 1).The mapper function does not modify the input data, but simply returns a new output list.

Figure 11 :
Figure11: Plots of the median errors in preserving the distances between pairs of data points versus the reduced dimension k in large scale microarray data.Here "RP" stands for the random projection and "PCA" stands for the principal component analysis.