Does the Ordinary Least Squares give the same estimation results as the LASSO with shrinkage parameter being zero? Does the Ordinary Least Squares give the same estimation results as the LASSO with shrinkage parameter being zero? If yes, is there any published papers or books to support this in theory?
 A: Short Answer
You don't need a published paper or book to support this. Setting $\lambda$ and thus the penalty to zero renders an equation that is identical to OLS: 
OLS:
$$\sum_{i=1}^n (y_i - \hat{y}_i)^2$$
LASSO
$$\sum_{i=1}^n (y_i - \hat{y}_i)^2 + \lambda\cdot\sum_{j=1}^p|\beta_j| $$
Setting $\lambda$ to zero leaves:
$$\sum_{i=1}^n (y_i - \hat{y}_i)^2 + 0\cdot\sum_{j=1}^p|\beta_j| = \sum_{i=1}^n (y_i - \hat{y}_i)^2$$
Which is OLS.

Elaborate Answer
OLS
In ordinary least squares, we estimate the coefficients ($\boldsymbol{\beta}$) of a $\color{darkblue}{\text{linear}}$ $\color{darkblue}{\text{model}}$:
$$\hat{y} = \color{darkblue}{\beta_0 + \sum_{j=i}^p \beta_j \cdot x_j}$$
by minimizing the squared $\color{darkred}{\text{difference}}$ $\color{darkred}{\text{from}}$ $\color{darkred}{\text{the}}$ $\color{darkred}{\text{observed values}}$. So if we have $p$ explanatory variables in: 
$$y_i = \color{darkblue}{\beta_0 + \sum_{j=i}^p \beta_j \cdot x_j} + \color{darkred}{\epsilon_i}$$
With $\beta_0$ the intercept and $\beta_j$ the coefficient for the $j$-th explanatory variable ($x_j$). We estimate this linear model such that the sum of these squared distances are minimized: 
$$\text{OLS} = \sum_{i=1}^n (\color{darkred}{y_i - \hat{y}_i})^2$$
Or visually (example using the circumference of a tree as a function of its age):

 

LASSO
The previous picture shows an example using just one explanatory variable. If we keep adding explanatory variables, eventually our model will end up overfitting (fitting the noise rather than the systematic component). Multiple explanatory variables are a little harder to visualize, so consider polynomial terms instead (this is just a mock example to illustrate overfitting; don't use LASSO to include many polynomial terms):  

Here we see that the sum of squared residuals is not a good minimization objective, because it will result in overfitting; the line follows the trend in the noise rather than just the systematic part of the process. So regularization methods (e.g. LASSO) address this problem by adding a $\color{darkgreen}{\text{penalty}}$ which constrains the coefficients (forces them to be small). LASSO in particular penalizes for the $\color{darkgreen}{\text{sum}}$ $\color{darkgreen}{\text{of}}$ $\color{darkgreen}{\text{absolute}}$ $\color{darkgreen}{\text{values}}$ $\color{darkgreen}{\text{of}}$ $\color{darkgreen}{\text{the}}$ $\color{darkgreen}{\text{coefficients}}$:  
$$\sum_{i=1}^n (\color{darkred}{y_i - \hat{y}_i})^2 + \lambda\cdot\color{darkgreen}{\sum_{j=1}^p|\beta_j|} $$
Note that this is just an addition of a penalty to OLS. By setting $\lambda$ to zero, the penalty becomes zero and what remains is just the original objective (OLS):
$$\sum_{i=1}^n (\color{darkred}{y_i - \hat{y}_i})^2 + 0\cdot\color{darkgreen}{\sum_{j=1}^p|\beta_j|} = \sum_{i=1}^n (\color{darkred}{y_i - \hat{y}_i})^2$$
As for your question about optimization in the comments, as SmallChess pointed out: both OLS and LASSO can be estimated with gradient descent. Only OLS can be estimated with maximum likelihood. If you want to be sure they yield the same results at $\lambda = 0$, try solving OLS with gradient descent.
