Why is the L2 regularization equivalent to Gaussian prior? I keep reading this and intuitively I can see this but how does one go from L2 regularization to saying that this is a Gaussian Prior analytically? Same goes for saying L1 is equivalent to a Laplacean prior. 
Any further references would be great. 
 A: Let us imagine that you want to infer some parameter $\beta$ from some observed input-output pairs $(x_1,y_1)\dots,(x_N,y_N)$. Let us assume that the outputs are linearly related to the inputs via $\beta$ and that the data are corrupted by some noise $\epsilon$:
$$y_n = \beta x_n + \epsilon,$$
where $\epsilon$ is Gaussian noise with mean $0$ and variance $\sigma^2$.
This gives rise to a Gaussian likelihood:
$$\prod_{n=1}^N \mathcal{N}(y_n|\beta x_n,\sigma^2).$$
Let us regularise parameter $\beta$ by imposing the Gaussian prior $\mathcal{N}(\beta|0,\lambda^{-1}),$ where $\lambda$ is a strictly positive scalar ($\lambda$ quantifies of by how much we believe that $\beta$ should be close to zero, i.e. it controls the strength of the regularisation).
Hence, combining the likelihood and the prior we simply have:
$$\prod_{n=1}^N \mathcal{N}(y_n|\beta x_n,\sigma^2) \mathcal{N}(\beta|0,\lambda^{-1}).$$
Let us take the logarithm of the above expression. Dropping some constants we get:
$$\sum_{n=1}^N -\frac{1}{\sigma^2}(y_n-\beta x_n)^2 - \lambda \beta^2 + \mbox{const}.$$
If we maximise the above expression with respect to $\beta$, we get the so called maximum a-posteriori estimate for $\beta$, or MAP estimate for short. In this expression it becomes apparent why the Gaussian prior can be interpreted as a L2 regularisation term.

The relationship between the L1 norm and the Laplace prior can be understood in the same fashion. Instead of a Gaussian prior, multiply your likelihood with a Laplace prior and then take the logarithm.

A good reference (perhaps slightly advanced) detailing both issues is the paper "Adaptive Sparseness for Supervised Learning", which currently does not seem easy to find online. Alternatively look at "Adaptive Sparseness using Jeffreys Prior". Another good reference is "On Bayesian classification with Laplace priors".
A: To put the equivalence more precisely:
Optimizing model weights to minimize a squared error loss function with L2 regularization is equivalent to finding the weights that are most likely under a posterior distribution evaluated using Bayes rule, with a zero-mean independent Gaussian weights prior
Proof:
The loss function as described above would be given by
$$
L = \underbrace{\Big[ \sum_{n=1}^{N} (y^{(n)} - f_{\mathbf{w}}(\mathbf{x}^{(n)}))^{2} \Big] }_{Original \; loss \; function} +  \underbrace{\lambda \sum_{i=1}^{K} w_{i}^{2}}_{L_{2} \; loss}
$$
Note that the distribution for a multivariate Gaussian is
$$
\mathcal{N}(\mathbf{x}; \mathbf{\mu}, \Sigma) = \frac{1}{(2 \pi)^{D/2}|\Sigma|^{1/2}} \exp\Big(-\frac{1}{2} (\mathbf{x} -\mathbf{\mu})^{\top} \Sigma^{-1} (\mathbf{x} -\mathbf{\mu})\Big)
$$
Using Bayes rule, we have that
$$
\begin{split}
p(\mathbf{w}|\mathcal{D}) &= \frac{p(\mathcal{D}|\mathbf{w}) \; p(\mathbf{w})}{p(\mathcal{D})}\newline
 &\propto p(\mathcal{D}|\mathbf{w}) \; p(\mathbf{w})\newline
&\propto \Big[ \prod_{n}^{N} \mathcal{N}(y^{(n)}; f_{\mathbf{w}}(\mathbf{x}^{(n)}), \sigma_{y}^{2})\Big] \; \mathcal{N}(\mathbf{w}; \mathbf{0}, \sigma_{\mathbf{w}}^{2} \mathbb{I})\newline
&\propto \prod_{n}^{N} \mathcal{N}(y^{(n)};f_{\mathbf{w}}(\mathbf{x}^{(n)}) , \sigma_{y}^{2}) \prod_{i=1}^{K} \mathcal{N}(w_{i}; \, 0, \, \sigma_{\mathbf{w}}^{2}) \newline
\end{split}
$$
Where we are able to split the multi-dimensional Guassian into a product, because the covariance is a multiple of the identity matrix. 
Take negative log probability
$$
\begin{split}
-\log \big[p(\mathbf{w}|\mathcal{D}) \big] &= -\sum_{n=1}^{N} \log \big[\mathcal{N}(y^{(n)}; f_{\mathbf{w}}(\mathbf{x}^{(n)}), \sigma_{y}^{2}) \big] - \sum_{i=1}^{K} \log \big[ \mathcal{N}(w_{i}; \, 0, \, \sigma_{\mathbf{w}}^{2}) \big] + const. \newline
&= \frac{1}{2\sigma_{y}^{2}} \sum_{n=1}^{N} \big(y^{(n)} - f_{\mathbf{w}}(\mathbf{x}^{(n)})\big)^{2} + \frac{1}{2\sigma_{\mathbf{w}}^{2}} \sum_{i=1}^{K} w_{i}^{2}   + const. \newline
\end{split}
$$
We can of course drop the constant, and multiply by any amount without fundamentally affecting the loss function. (constant does nothing,  multiplication effectively scales the learning rate. Will not affect the location of minima)So we can see that the negative log probability of the posterior distribution is an equivalent loss function to the L2 regularized square error loss function. 
This equivelance is general and holds for any parameterized function of weights - not just linear regression as seems to be implied above. 
A: There are two characteristics of Bayesian modeling that need to be emphasized, when discussing the equivalance of certain penalized maximum likelihood estimation and Bayesian procedures.


*

*In the Bayesian framework, the prior is selected based on specifics of the problem and is not motivated by computational expediency.  Hence Bayesians use a variety of priors including the now-popular horseshoe prior for sparse predictor problems, and don't need to rely so much on priors that are equivalent to L1 or L2 penalties.

*With a full Bayesian approach you have access to all inferential procedures when you're done.  For example you can quantify evidence for large regression coefficients and you can get credible intervals on regression coefficients and overall predicted values.  In the frequentist framework, once you choose penalization you lose all of the inferential machine.

A: First notice that median minimizes the L1 norm (see here or here for learning more on L1 and L2)
$$ \DeclareMathOperator*{\argmin}{arg\,min}
\text{median}(x) = \argmin_s \sum_i |x_i - s|^1  $$
while mean minimizes L2
$$ \text{mean}(x) = \argmin_s \sum_i |x_i - s|^2  $$
now, recall that Normal distributions' $\mu$ parameter can be estimated using sample mean, while the MLE estimator for Laplace distribution $\mu$ parameter is median. So using Normal distribution is equivalent to L2 norm optimization and using Laplace distribution, to using L1 optimization. In practice you can think of it as that median is less sensitive to outliers than mean, and the same, using fatter-tailed Laplace distribution as a prior makes your model less prone to outliers, than using Normal distribution.

Hurley, W. J. (2009) An Inductive Approach to Calculate the MLE for the Double Exponential Distribution. Journal of Modern Applied Statistical Methods: 8(2), Article 25.
A: For a linear model with multivariate normal prior and multivariate normal likelihood, you end up with a multivariate normal posterior distribution in which the mean of the posterior (and maximum a posteriori model) is exactly what you would obtain using Tikhonov regularized ($L_{2}$ regularized) least squares with an appropriate regularization parameter.  
Note that there is a more fundamental difference in that the Bayesian posterior is a probability distribution, while the Tikhonov regularized least squares solution is a specific point estimate.  
This is discussed in many textbooks on Bayesian methods for inverse problems, See for example:
http://www.amazon.com/Inverse-Problem-Methods-Parameter-Estimation/dp/0898715725/
http://www.amazon.com/Parameter-Estimation-Inverse-Problems-Second/dp/0123850487/
Similarly, if you have a Laplacian prior and a multivariate normal likelihood, then the maximum of the posterior distribution occurs at a point that you could get by solving an $L_{1}$ regularized least squares problem.  
A: For a regression problem with $k$ variables (w/o intercept) you do OLS as
$$\min_{\beta} (y - X \beta)' (y - X \beta)$$
In regularized regression with $L^p$ penalty you do
$$\min_{\beta} (y - X \beta)' (y - X \beta) + \lambda \sum_{i=1}^k |\beta_i|^p $$
We can equivalently do (note the sign changes)
$$\max_{\beta} -(y - X \beta)' (y - X \beta) - \lambda \sum_{i=1}^k |\beta_i|^p $$
This directly relates to the Bayesian principle of
$$posterior \propto likelihood \times prior$$
or equivalently (under regularity conditions)
$$log(posterior) \sim log(likelihood) + log(penalty)$$
Now it is not hard to see which exponential family distribution corresponds to which penalty type.
