On the importance of the i.i.d. assumption in statistical learning In statistical learning, implicitly or explicitly, one always assumes that the training set $\mathcal{D} = \{ \bf {X}, \bf{y} \}$ is composed of $N$ input/response tuples $({\bf{X}}_i,y_i)$ that are independently drawn from the same joint distribution $\mathbb{P}({\bf{X}},y)$ with
$$ p({\bf{X}},y) = p( y \vert {\bf{X}}) p({\bf{X}}) $$
and $p( y \vert {\bf{X}})$ the relationship we are trying to capture through a particular learning algorithm. Mathematically, this i.i.d. assumption writes:
\begin{gather}
({\bf{X}}_i,y_i) \sim \mathbb{P}({\bf{X}},y), \forall i=1,...,N \\
({\bf{X}}_i,y_i) \text{ independent of } ({\bf{X}}_j,y_j), \forall i \ne j \in \{1,...,N\}
\end{gather}
I think we can all agree that this assumption is rarely satisfied in practice, see this related SE question and the wise comments of @Glen_b and @Luca.
My question is therefore: 

Where exactly does the i.i.d. assumption becomes critical in practice? 

[Context] 
I'm asking this because I can think of many situations where such a stringent assumption is not needed to train a certain model (e.g. linear regression methods), or at least one can work around the i.i.d. assumption and obtain robust results. Actually the results will usually stay the same, it is rather the inferences that one can draw that will change (e.g. heteroskedasticity and autocorrelation consistent HAC estimators in linear regression: the idea is to re-use the good old OLS regression weights but to adapt the finite-sample behaviour of the OLS estimator to account for the violation of the Gauss-Markov assumptions). 
My guess is therefore that the i.i.d. assumption is required not to be able to train a particular learning algorithm, but rather to guarantee that techniques such as cross-validation can indeed be used to infer a reliable measure of the model's capability of generalising well, which is the only thing we are interested in at the end of the day in statistical learning because it shows that we can indeed learn from the data. Intuitively, I can indeed understand that using cross-validation on dependent data could be optimistically biased (as illustrated/explained in this interesting example). 
For me i.i.d. has thus nothing to do with training a particular model but everything to do with that model's generalisability. This seems to agree with a paper I found by Huan Xu et al, see "Robustness and Generalizability for Markovian Samples" here. 
Would you agree with that?
[Example] 
If this can help the discussion, consider the problem of using the LASSO algorithm to perform a smart selection amongst $P$ features given $N$ training samples $({\bf{X}}_i,y_i)$ with $\forall i=1,...,N$
$$ {\bf{X}}_i=[X_{i1},...,X_{iP}] $$
We can further assume that:


*

*The inputs ${\bf{X}}_i$ are dependent hence leading to a violation of the i.i.d. assumption (e.g. for each feature $j=1,..,P$ we observe a $N$ point time series, hence introducing temporal auto-correlation)

*The conditional responses $y_i \vert {\bf{X}}_i$ are independent.

*We have $P \gg N$. 


In what way(s) does the violation of the i.i.d. assumption can pose problem in that case assuming we plan to determine the LASSO penalisation coefficient $\lambda$ using a cross-validation approach (on the full data set) + use a nested cross-validation to get a feel for the generalisation error of this learning strategy (we can leave the discussion concerning the inherent pros/cons of the LASSO aside, except if it is useful).
 A: The i.i.d. assumption about the pairs $(\mathbf{X}_i, y_i)$, $i = 1, \ldots, N$, is often made in statistics and in machine learning. Sometimes for a good reason, sometimes out of convenience and sometimes just because we usually make this assumption. To satisfactorily answer if the assumption is really necessary, and what the consequences are of not making this assumption, I would easily end up writing a book (if you ever easily end up doing something like that). Here I will try to give a brief overview of what I find to be the most important aspects.
A fundamental assumption
Let's assume that we want to learn a probability model of $y$ given $\mathbf{X}$, which we call $p(y \mid \mathbf{X})$. We do not make any assumptions about this model a priory, but we will make the minimal assumption that such a model exists such that

*

*the conditional distribution of $y_i$ given $\mathbf{X}_i$ is $p(y_i \mid \mathbf{X}_i)$.

What is worth noting about this assumption is that the conditional distribution of $y_i$ depends on $i$ only through $\mathbf{X}_i$. This is what makes the model useful, e.g. for prediction. The assumption holds as a consequence of the identically distributed part under the i.i.d. assumption, but it is weaker because we don't make any assumptions about the $\mathbf{X}_i$'s.
In the following the focus will mostly be on the role of independence.
Modelling
There are two major approaches to learning a model of $y$ given $\mathbf{X}$. One approach is known as discriminative modelling and the other as generative modelling.

*

*Discriminative modelling: We model $p(y \mid \mathbf{X})$ directly, e.g. a logistic regression model, a neural network, a tree or a random forest. The working modelling assumption will typically be that the $y_i$'s are conditionally independent given the $\mathbf{X}_i$'s, though estimation techniques relying on subsampling or bootstrapping make most sense under the i.i.d. or the weaker exchangeability assumption (see below). But generally, for discriminative modelling we don't need to make distributional assumptions about the $\mathbf{X}_i$'s.

*Generative modelling: We model the joint distribution, $p(\mathbf{X}, y)$, of $(\mathbf{X}, y)$ typically by modelling the conditional distribution $p(\mathbf{X} \mid y)$ and the marginal distribution $p(y)$. Then we use Bayes's formula for computing $p(y \mid \mathbf{X})$. Linear discriminant analysis and naive Bayes methods are examples. The working modelling assumption will typically be the i.i.d. assumption.

For both modelling approaches the working modelling assumption is used to derive or propose learning methods (or estimators). That could be by maximising the (penalised) log-likelihood, minimising the empirical risk or by using Bayesian methods. Even if the working modelling assumption is wrong, the resulting method can still provide a sensible fit of $p(y \mid \mathbf{X})$.
Some techniques used together with discriminative modelling, such as bagging (bootstrap aggregation), work by fitting many models to data sampled randomly from the dataset. Without the i.i.d. assumption (or exchangeability) the resampled datasets will not have a joint distribution similar to that of the original dataset. Any dependence structure has become "messed up" by the resampling. I have not thought deeply about this, but I don't see why that should necessarily break the method as a method for learning $p(y \mid \mathbf{X})$. At least not for methods based on the working independence assumptions. I am happy to be proved wrong here.
Consistency and error bounds
A central question for all learning methods is whether they result in models close to $p(y \mid \mathbf{X})$. There is a vast theoretical literature in statistics and machine learning dealing with consistency and error bounds. A main goal of this literature is to prove that the learned model is close to $p(y \mid \mathbf{X})$ when $N$ is large. Consistency is a qualitative assurance, while error bounds provide (semi-) explicit quantitative control of the closeness and give rates of convergence.
The theoretical results all rely on assumptions about the joint distribution of the observations in the dataset. Often the working modelling assumptions mentioned above are made (that is, conditional independence for discriminative modelling and i.i.d. for generative modelling). For discriminative modelling, consistency and error bounds will require that the $\mathbf{X}_i$'s fulfil certain conditions. In classical regression one such condition is that $\frac{1}{N} \mathbb{X}^T \mathbb{X} \to \Sigma$ for $N \to \infty$, where $\mathbb{X}$ denotes the design matrix with rows $\mathbf{X}_i^T$. Weaker conditions may be enough for consistency. In sparse learning another such condition is the restricted eigenvalue condition, see e.g. On the conditions used to prove oracle results for the Lasso. The i.i.d. assumption together with some technical distributional assumptions imply that some such sufficient conditions are fulfilled with large probability, and thus the i.i.d. assumption may prove to be a sufficient but not a necessary assumption to get consistency and error bounds for discriminative modelling.
The working modelling assumption of independence may be wrong for either of the modelling approaches. As a rough rule-of-thumb one can still expect consistency if the data comes from an ergodic process, and one can still expect some error bounds if the process is sufficiently fast mixing. A precise mathematical definition of these concepts would take us too far away from the main question. It is enough to note that there exist dependence structures besides the i.i.d. assumption for which the learning methods can be proved to work as $N$ tends to infinity.
If we have more detailed knowledge about the dependence structure, we may choose to replace the working independence assumption used for modelling with a model that captures the dependence structure as well. This is often done for time series. A better working model may result in a more efficient method.
Model assessment
Rather than proving that the learning method gives a model close to $p(y \mid \mathbf{X})$ it is of great practical value to obtain a (relative) assessment of "how good a learned model is". Such assessment scores are comparable for two or more learned models, but they will not provide an absolute assessment of how close a learned model is to $p(y \mid \mathbf{X})$. Estimates of assessment scores are typically computed empirically based on splitting the dataset into a training and a test dataset or by using cross-validation.
As with bagging, a random splitting of the dataset will "mess up" any dependence structure. However, for methods based on the working independence assumptions, ergodicity assumptions weaker than i.i.d. should be sufficient for the assessment estimates to be reasonable, though standard errors on these estimates will be very difficult to come up with.
[Edit: Dependence among the variables will result in a distribution of the learned model that differs from the distribution under the i.i.d. assumption. The estimate produced by cross-validation is not obviously related to the generalization error. If the dependence is strong, it will most likely be a poor estimate.]
Summary (tl;dr)
All the above is under the assumption that there is a fixed conditional probability model, $p(y \mid \mathbf{X})$. Thus there cannot be trends or sudden changes in the conditional distribution not captured by $\mathbf{X}$.
When learning a model of $y$ given $\mathbf{X}$, independence plays a role as

*

*a useful working modelling assumption that allows us to derive learning methods

*a sufficient but not necessary assumption for proving consistency and providing error bounds

*a sufficient but not necessary assumption for using random data splitting techniques such as bagging for learning and cross-validation for assessment.

To understand precisely what alternatives to i.i.d. that are also sufficient is non-trivial and to some extent a research subject.
A: In my opinion there are two rather mundane reasons why the i.i.d. assumption is important in statistical learning (or statistics in general).


*

*Lots of behind the scenes mathematics depend on this assumption. If you want to prove that your learning method actually works for more than one data set, i.i.d. assumption will crop up eventually. It is possible to avoid it, but mathematics becomes several times harder. 

*If you want to learn something from data, you need to assume that there is something to learn. Learning is impossible if every data point is generated by different mechanism. So it is essential to assume that something unifies given data set. If we assume that data is random, then this something is naturally a probability distribution, because probability distribution encompasses all information about the random variable. 
So if we have data $x_1,...,x_n$ ($x_i$ can be either a vector or scalar), we assume that it comes from distribution $F_n$: 
$$(x_1,...,x_n)\sim F_n.$$ 
Here we have a problem. We need to ensure that $F_n$ is related to $F_m$, for different $n$ and $m$, otherwise we have the initial problem, that every data point is generated differently. The second problem is that although we have $n$ data points, we basically have one data point for estimating $F_n$, because $F_n$ is $n$-variate probability distribution. The most simple solution for these two problems is an i.i.d assumption. With it $F_n=F^n,$ where $x_i\sim F$. We get very clear relationship between $F_n$ and $F_m$ and we have $n$ data points to estimate one $F$. There are other ways these two problems are solved, but it is essential to note that every statistical learning method needs to solve this problem and it so happens that i.i.d. assumption is by far the most uncomplicated way to do it.
A: The only place where one can safely ignored iid is in undergraduate statistics and machine learning courses. You have written that:

one can work around the i.i.d. assumption and obtain robust results. Actually the results will usually stay the same, it is rather the inferences that one can draw that will change...

This is only true if the functional form of the models is assumed to be basically correct. But, such an assumption is even less plausible than iid.
There are at least two ways in which iid is critically important in terms of applied modeling:


*

*It is an explicit assumption in most statistical inference, as you note in your question. In most real-world modeling, at some stage we need to use inference to test the specification, such as during variable selection and model comparison. So, while each particular model fit may be OK despite iid violations, you can end up choosing the wrong model anyway.

*I find that thinking through violations of iid is a useful way to think about the data generating mechanism, which in turn helps me think about the appropriate specification of a model a priori. Two examples:


*

*If the data is clustered, this is a violation of iid. A remedy to this may be a mixture model. The inference I will draw from a mixture models is generally completely different to that which I draw from OLS.

*Non-linear relationships between the dependent and independent variables often show up when inspecting residuals as a part of investigating iid.



Of course, in pretty much ever model that I have ever built, I have failed in my quest to reduce the distribution of the residuals to anything close to a truly normal distribution. But, nevertheless, I always gain a lot by trying really, really, hard to do it.
A: What i.i.d. assumption states is that random variables are independent and identically distributed. You can formally define what does it mean, but informally it says that all the variables provide the same kind of information independently of each other (you can read also about related exchangeability).
From the abstract ideas let's jump for a moment to concrete example: in most cases your data can be stored in a matrix, with observations row-wise and variables column-wise. If you assume your data to be i.i.d., then it means for you that you need to bother only about relations between columns and do not have to bother about relations between rows. If you bothered about both then you would model dependence of columns on columns and rows on rows, i.e. everything on everything. It is very hard to make simplifications and build a statistical model of everything depending on everything.
You correctly noticed that exchengeability makes it possible for us to use methods such as cross-validation, or bootstrap, but it also makes it possible to use central limit theorem and it enables us to make simplifications helpful for modeling (thinking in column-wise terms).
As you noticed in the LASSO example, independence assumption is often softened to conditional independence. Even in such case we need independent and identically distributed "parts". Similar, softer assumption is often made for time-series models, that you mentioned, that assume stationarity (so there is dependence but there is also a common distribution and series stabilizes over time -- again "i.i.d." parts). It is a matter of observing a number of similar things that carry the same idea about some general phenomenon. If we have a number of distinct and dependent things we cannot make any generalizations.
What you have to remember is that this is only an assumption, we are not strict about it. It is about having enough things that all, independently, convey similar information about some common phenomenon. If the things influenced each other they would obviously convey similar information so they wouldn't be that useful.
Imagine that you wanted to learn about abilities of children in a classroom, so you give them some tests. You could use the test results as an indicator of the abilities of kids only if they did them by themselves, independently of each other. If they interacted then you'd probably measure abilities of the most clever kid, or the most influential one. It does not mean that you need to assume that there was no interaction, or dependence, between kids whatsoever, but simply that they did the tests by themselves. The kids also need to be "identically distributed", so they cannot come from different countries, speak different languages, be in different ages since it will make it hard to interpret the results (maybe they did not understand the questions and answered randomly). If you can assume that your data is i.i.d. then you can focus on building a general model. You can deal with non-i.i.d. data but then you have to worry about "noise" in your data much more.

Besides your main question you are also asking about cross-validation with non-i.i.d. data. While you seem to understate the importance of i.i.d. assumption, at the same time you overstate the problems of not meeting this assumption poses for cross-validation. There are multiple ways how we can deal with such data when using resampling methods like bootstrap, or cross-validation. If you are dealing with time-series you cannot assume that the values are independent, so taking the random fraction of values would be a bad idea because it would ignore the autocorrelated structure of the data. Because of that, with time-series we commonly use one step ahead cross-validation, i.e. you take part of the series to predict next value (not used for modeling). Similarly, if your data has clustered structure, you sample whole clusters to preserve the nature of the data. So as with modeling, we can deal with non-i.i.d.-sness also when doing cross-validation, but we need to adapt our methods to the nature of the data since methods designed for i.i.d. data do not apply in such cases.
A: I would like to stress that in some circumstances, the data are not i.i.d. and statistical learning is still possible. It is crucial to have an identifiable model for the joint distribution of all observations; if the observations are i.i.d. then this joint distribution is easily obtained from the marginal distribution of single observations. But in some cases, the joint distribution is given directly, without resorting to a marginal distribution.
A widely used model in which the observations are not i.i.d. is the linear mixed model:
$$\let\epsilon\varepsilon Y = X \alpha + Z u + \epsilon $$
with $\def\R{\mathbb{R}}Y \in \R^n$, $X \in \R^{n\times p}$, $\alpha \in \R^p$, $Z \in \R^{n\times q}$, $u \in \R^q$, and $\epsilon\in\R^n$. The (design) matrix $X$ and $Z$ are considered as fixed, $\alpha$ is a vector of parameters, $u$ is a random vector $\def\N{\mathcal{N}} u\sim \N(0,\tau I_q)$ and $\epsilon \sim \N(0,\sigma^2 I_n)$, $\tau$ and $\sigma^2$ being parameters of the model.
This model is best expressed by giving the distribution of $Y$:
$$Y \sim \N(X\alpha, \tau ZZ' + \sigma^2 I_n).$$
The parameters to be learned are $\alpha$, $\tau$, $\sigma^2$. A single vector $Y$ of dimension $n$ is observed; its components are not i.i.d. 
A: One area where i.i.d. assumption is critical in practice (other that inference) is data collection. If you do not collect data in a random manner then you will have a sampling bias and your data will not be a good representation of the underlying model.
