What is the difference between linear regression on y with x and x with y? The Pearson correlation coefficient of x and y is the same, whether you compute pearson(x, y) or pearson(y, x). This suggests that doing a linear regression of y given x or x given y should be the same, but I don't think that's the case. 
Can someone shed light on when the relationship is not symmetric, and how that relates to the Pearson correlation coefficient (which I always think of as summarizing the best fit line)?
 A: On questions like this it's easy to get caught up on the technical issues, so I'd like to focus specifically on the question in the title of the thread which asks: What is the difference between linear regression on y with x and x with y?
Consider for a moment a (simplified) econometric model from human capital theory (the link goes to an article by Nobel Laureate Gary Becker). Let's say we specify a model of the following form:
\begin{equation} \text{wages} = b_{0} + b_{1}~\text{years of education}  + \text{error}
\end{equation}
This model can be interpreted as a causal relationship between wages and education. Importantly, causality in this context means the direction of causality runs from education to wages and not the other way round. This is implicit in the way the model has been formulated; the dependent variable is wages and the independent variable is years of education.
Now, if we make a reversal of the econometric equation (that is, change y on x to x on y), such that the model becomes
\begin{equation} \text{years of education} = b_{0} + b_{1}~\text{wages} + \text{error}
\end{equation}
then implicit in the formulation of the econometric equation is that we are saying that the direction of causality runs from wages to education.
I'm sure you can think of more examples like this one (outside the realm of economics too), but as you can see, the interpretation of the model can change quite significantly when we switch from regressing y on x to x on y.
So, to the answer the question: What is the difference between linear regression on y with x and x with y?, we can say that the interpretation of the regression equation changes when we regress x on y instead of y on x. We shouldn't overlook this point because a model that has a sound interpretation can quickly turn into one which makes little or no sense.
A: Expanding on @gung's excellent answer:  
In a simple linear regression the absolute value of Pearson's $r$ can be seen as the geometric mean of the two slopes we obtain if  we regress $y$ on $x$ and $x$ on $y$, respectively:
$$\sqrt{{\hat{\beta}_1}_{y\,on\,x} \cdot {\hat{\beta}_1}_{x\,on\,y}} = \sqrt{\frac{\text{Cov}(x,y)}{\text{Var}(x)} \cdot \frac{\text{Cov}(y,x)}{\text{Var}(y)}} = \frac{|\text{Cov}(x,y)|}{\text{SD}(x) \cdot \text{SD}(y)} = |r|
$$
We can obtain $r$ directly using
$$r = sign({\hat{\beta}_1}_{y\,on\,x}) \cdot \sqrt{{\hat{\beta}_1}_{y\,on\,x} \cdot {\hat{\beta}_1}_{x\,on\,y}}
$$ 
or 
$$r = sign({\hat{\beta}_1}_{x\,on\,y}) \cdot \sqrt{{\hat{\beta}_1}_{y\,on\,x} \cdot {\hat{\beta}_1}_{x\,on\,y}}
$$
Interestingly, by the AM–GM inequality, it follows that the absolute value of the arithmetic mean of the two slope coefficients is greater than (or equal to) the absolute value of Pearson's $r$:
$$
|\frac{1}{2} \cdot ({\hat{\beta}_1}_{y\,on\,x} + {\hat{\beta}_1}_{x\,on\,y})| \geq \sqrt{{\hat{\beta}_1}_{y\,on\,x} \cdot {\hat{\beta}_1}_{x\,on\,y}} = |r|
$$
A: There is a very interesting phenomenon about this topic. After exchanging x and y, although the regression coefficient changes, but the t-statistic/F-statistic and significance level for the coefficient don't change. This is also true even in  multiple regression, where we exchange y with one of the independent variables.
It is due to a delicate relation between the F-statistic and (partial) correlation coefficient. That relation really touches the core of linear model theory.There are more details about this conclusion in my notebook: Why exchange y and x has no effect on p
A: The relation is not symmetric because we are solving two different optimisation problems. $\textbf{ Doing regression of $y$ given $x$}$ can be written as solving the following problem:
$$\min_b \mathbb E(Y - bX)^2$$
whereas for $\textbf{doing regression of $x$ given $y$}$:
$$\min_b \mathbb E(X - bY)^2$$, which can be rewritten as:
$$\min_b \frac{1}{b^2} \mathbb E(Y - bX)^2$$
It is also important to note that, two different-looking problems may have the same solution.
A: The best way to think about this is to imagine a scatterplot of points with $y$ on the vertical axis and $x$ represented by the horizontal axis.  Given this framework, you see a cloud of points, which may be vaguely circular, or may be elongated into an ellipse.  What you are trying to do in regression is find what might be called the 'line of best fit'.  However, while this seems straightforward, we need to figure out what we mean by 'best', and that means we must define what it would be for a line to be good, or for one line to be better than another, etc.  Specifically, we must stipulate a loss function.  A loss function gives us a way to say how 'bad' something is, and thus, when we minimize that, we make our line as 'good' as possible, or find the 'best' line.  
Traditionally, when we conduct a regression analysis, we find estimates of the slope and intercept so as to minimize the sum of squared errors.  These are defined as follows: 
$$
SSE=\sum_{i=1}^N(y_i-(\hat\beta_0+\hat\beta_1x_i))^2
$$  
In terms of our scatterplot, this means we are minimizing the (sum of the squared) vertical distances between the observed data points and the line.  

On the other hand, it is perfectly reasonable to regress $x$ onto $y$, but in that case, we would put $x$ on the vertical axis, and so on.  If we kept our plot as is (with $x$ on the horizontal axis), regressing $x$ onto $y$ (again, using a slightly adapted version of the above equation with $x$ and $y$ switched) means that we would be minimizing the sum of the horizontal distances between the observed data points and the line.  This sounds very similar, but is not quite the same thing.  (The way to recognize this is to do it both ways, and then algebraically convert one set of parameter estimates into the terms of the other.  Comparing the first model with the rearranged version of the second model, it becomes easy to see that they are not the same.)  

Note that neither way would produce the same line we would intuitively draw if someone handed us a piece of graph paper with points plotted on it.  In that case, we would draw a line straight through the center, but minimizing the vertical distance yields a line that is slightly flatter (i.e., with a shallower slope), whereas minimizing the horizontal distance yields a line that is slightly steeper.
A correlation is symmetrical; $x$ is as correlated with $y$ as $y$ is with $x$.  The Pearson product-moment correlation can be understood within a regression context, however.  The correlation coefficient, $r$, is the slope of the regression line when both variables have been standardized first.  That is, you first subtracted off the mean from each observation, and then divided the differences by the standard deviation.  The cloud of data points will now be centered on the origin, and the slope would be the same whether you regressed $y$ onto $x$, or $x$ onto $y$ (but note the comment by @DilipSarwate below).  

Now, why does this matter?  Using our traditional loss function, we are saying that all of the error is in only one of the variables (viz., $y$).  That is, we are saying that $x$ is measured without error and constitutes the set of values we care about, but that $y$ has sampling error.  This is very different from saying the converse.  This was important in an interesting historical episode:  In the late 70's and early 80's in the US, the case was made that there was discrimination against women in the workplace, and this was backed up with regression analyses showing that women with equal backgrounds (e.g., qualifications, experience, etc.) were paid, on average, less than men.  Critics (or just people who were extra thorough) reasoned that if this was true, women who were paid equally with men would have to be more highly qualified, but when this was checked, it was found that although the results were 'significant' when assessed the one way, they were not 'significant' when checked the other way, which threw everyone involved into a tizzy.  See here for a famous paper that tried to clear the issue up.  

(Updated much later)  Here's another way to think about this that approaches the topic through the formulas instead of visually:  
The formula for the slope of a simple regression line is a consequence of the loss function that has been adopted.  If you are using the standard Ordinary Least Squares loss function (noted above), you can derive the formula for the slope that you see in every intro textbook.  This formula can be presented in various forms; one of which I call the 'intuitive' formula for the slope.  Consider this form for both the situation where you are regressing $y$ on $x$, and where you are regressing $x$ on $y$:
$$
\overbrace{\hat\beta_1=\frac{\text{Cov}(x,y)}{\text{Var}(x)}}^{y\text{ on } x}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\overbrace{\hat\beta_1=\frac{\text{Cov}(y,x)}{\text{Var}(y)}}^{x\text{ on }y}
$$
Now, I hope it's obvious that these would not be the same unless $\text{Var}(x)$ equals $\text{Var}(y)$.  If the variances are equal (e.g., because you standardized the variables first), then so are the standard deviations, and thus the variances would both also equal $\text{SD}(x)\text{SD}(y)$.  In this case, $\hat\beta_1$ would equal Pearson's $r$, which is the same either way by virtue of the principle of commutativity: 
$$
\overbrace{r=\frac{\text{Cov}(x,y)}{\text{SD}(x)\text{SD}(y)}}^{\text{correlating }x\text{ with }y}~~~~~~~~~~~~~~~~~~~~~~~~~~~\overbrace{r=\frac{\text{Cov}(y,x)}{\text{SD}(y)\text{SD}(x)}}^{\text{correlating }y\text{ with }x}
$$
A: This question can also be answered from a linear algebra perspective. Say you have a bunch of data points $(x,y)$. We want to find the line $y=mx+b$ that's closest to all our points (the regression line). 
As an example, say we have the points $(1,2),(2,4.5),(3,6),(4,7)$. We can look at this as a simultaneous equation problem:
\begin{align}
 & \underline{mx + b = y}\\
 & 1x + b = 2 \\ 
 & 2x + b = 4.5 \\
 & 3x + b = 6 \\ 
 & 4x + b = 7
\end{align}
In matrix form:
$$
    \left[\begin{matrix}
    1 & 1 \\
    2 & 1 \\
    3 & 1 \\
    4 & 1 
    \end{matrix}\right]    \left[\begin{matrix}
    x  \\
    b  \\
    \end{matrix}\right]=\left[\begin{matrix}
    2 \\
    4.5 \\
    6 \\
    7  
    \end{matrix}\right] 
$$
We see right away that $\vec{y}=(2,4.5,6,7)$ (the right hand side vector) is not in the span of the columns of our matrix, meaning we will not find an $(x,b)$ to solve our system.
The closest vector to $\vec{y}$ we can find in our column space is the projection $\vec p$ of $\vec{y}$ on the column space.
If we swap out $\vec{y}$ with its projection $\vec p$ on the column space, and solve our  system of equations for $\vec p$, we get the least squares solution, aka the regression line.
I.e. we can solve
$$
    \left[\begin{matrix}
    1 & 1 \\
    2 & 1 \\
    3 & 1 \\
    4 & 1 
    \end{matrix}\right]    \left[\begin{matrix}
    x  \\
    b  \\
    \end{matrix}\right]=\left[\begin{matrix}
    p_1 \\
    p_2 \\
    p_3 \\
    p_4  
    \end{matrix}\right] 
$$
to obtain the regression line $y=mx+b$ (here $m$ is the correlation coefficient normally called $\beta$).  
If you did $x=my+b$ instead, you'd have:
$$
    \left[\begin{matrix}
    2 & 1 \\
    4.5 & 1 \\
    6 & 1 \\
    7 & 1 
    \end{matrix}\right]    \left[\begin{matrix}
    y  \\
    b  \\
    \end{matrix}\right]=\left[\begin{matrix}
    1 \\
    2 \\
    3 \\
    4  
    \end{matrix}\right] 
$$
To find the regression line, we'd have to solve this system using the projection $\vec r$ of $\vec x = (1,2,3,4)$ on to the column space of our new matrix. 
That is, we swap $(1,2,3,4)$ with its projection $(r_1,r_2,r_3,r_4)$ on the span of $(2,4.5,6,7)$ and $(1,1,1,1)$ and solve the system. You can solve it by hand if you want to and compare it to a least squares solution found by a computer.
The idea that the regression of y given x or x given y should be the same, is equivalent to asking if $\vec p=\vec r$ in linear algebra terms.
We know that $\vec p$ is in $span (\vec x,\vec b)$ and $\vec r$ is in $span (\vec y,\vec b)$. We known that $\vec x \neq c \vec y$ since this is what motivated us to look for a regression line in the first place.
Therefore, the intersection of $span (\vec x,\vec b)$ and $span (\vec y,\vec b)$ is $c \vec b$.
So if $\vec p=\vec r$, then $\vec p=\vec r = c \vec b$.
What type of line is  $c\vec b = c(1,1,1,\dots)$? On the plane, it's $y=x$. It's the line that goes out 45° from the axes of your plot.
Most of the time our regression lines will not be of the $y=x$ type. So we can see how regression is usually not symmetric.
The correlation is symmetric however. From a linear algebra perspective the correlation (aka pearson(x,y)) is  $\cos(\theta)$ where $\theta$ is the angle between $\vec x$ and $\vec y$. 
In the example, the correlation/pearson(x,y) is the $\cos(\theta)$ of $(1,2,3,4)$ and $(2,4.5,6,7)$.
Clearly the angle between  $\vec x$ and $\vec y$ is equal to the angle between $\vec y$ and $\vec x$, so the correlation must be too.
A: I'm going to illustrate the answer with some R code and output.
First, we construct a random normal distribution, y, with a mean of 5 and a SD of 1:
y <- rnorm(1000, mean=5, sd=1)

Next, I purposely create a second random normal distribution, x, which is simply 5x the value of y for each y:
x <- y*5

By design, we have perfect correlation of x and y:
cor(x,y)
[1] 1
cor(y,x)
[1] 1

However, when we do a regression, we are looking for a function that relates x and y so the results of the regression coefficients depend on which one we use as the dependent variable, and which we use as the independent variable. In this case, we don't fit an intercept because we made x a function of y with no random variation:
lm(y~x-1)
Call:
lm(formula = y ~ x - 1)

Coefficients:
  x  
0.2

lm(x ~ y-1)
Call:
lm(formula = x ~ y - 1)

Coefficients:
y  
5  

So the regressions tell us that y=0.2x and that x=5y, which of course are equivalent. The correlation coefficient is simply showing us that there is an exact match in unit change levels between x and y, so that (for example) a 1-unit increase in y always produces a 0.2-unit increase in x.
A: The insight that since Pearson's correlation is the same whether we do a regression of x against y, or y against x is a good one, we should get the same linear regression is a good one.  It is only slightly incorrect, and we can use it to understand what is actually occurring.
This is the equation for a line, which is what we are trying to get from our regression

The equation for the slope of that line is driven by Pearson's correlation

This is the equation for Pearson's correlation.  It is the same whether we are regressing x against y or y against x

However when we look back at our second equation for slope, we see that Pearson's correlation is not the only term in that equation.  If we are calculating y against x, we also have the sample standard deviation of y divided by the sample standard deviation of x.  If we were to calculate the regression of x against y we would need to invert those two terms.
A: Well, it's true that for a simple bivariate regression, the linear correlation coefficient and R-square will be the same for both equations.  But the slopes will be $rS_y/S_x$ or $rS_x/S_y$ , which are not reciprocals of each other, unless $r = 1$.
