[proofplan]
We compare $f$ not directly with its interpolant, but with an arbitrary polynomial $p \in \mathcal{P}_N$. Since interpolation fixes every polynomial of degree at most $N$, the error decomposes as the approximation error $f-p$ plus the interpolation of that same error. The operator norm of $I_N$ on $C([-1,1])$ is bounded by the Lebesgue constant because the Lagrange formula expresses $I_N g$ as a signed weighted sum of the nodal values of $g$. Taking the infimum over all $p \in \mathcal{P}_N$ gives the desired inequality.
[/proofplan]
[step:Express the interpolation operator through the Lagrange basis]
For every $g \in C([-1,1])$, the interpolant $I_N g \in \mathcal{P}_N$ is given by
\begin{align*}
(I_N g)(x)=\sum_{j=0}^N g(x_j)\ell_j(x)
\end{align*}
for every $x \in [-1,1]$.
Indeed, the right-hand side is a polynomial in $\mathcal{P}_N$, and at each interpolation node $x_i$ it has value
\begin{align*}
\sum_{j=0}^N g(x_j)\ell_j(x_i)=\sum_{j=0}^N g(x_j)\delta_{ij}=g(x_i).
\end{align*}
The uniqueness of polynomial interpolation at distinct nodes then identifies this polynomial with $I_N g$.
[/step]
[step:Bound the interpolation operator by the Lebesgue constant]
Let $\|\cdot\|_\infty$ denote the supremum norm on $C([-1,1])$. For every $g \in C([-1,1])$ and every $x \in [-1,1]$, the Lagrange formula and the triangle inequality give
\begin{align*}
|(I_N g)(x)|\leq \sum_{j=0}^N |g(x_j)| |\ell_j(x)|.
\end{align*}
Since $|g(x_j)| \leq \|g\|_\infty$ for each $j$, this implies
\begin{align*}
|(I_N g)(x)|\leq \|g\|_\infty \sum_{j=0}^N |\ell_j(x)|.
\end{align*}
Taking the maximum over $x \in [-1,1]$ gives
\begin{align*}
\|I_N g\|_\infty \leq \lambda_N \|g\|_\infty.
\end{align*}
Thus the operator norm of $I_N:C([-1,1])\to C([-1,1])$ satisfies $\|I_N\| \leq \lambda_N$.
[guided]
The point of introducing the Lebesgue constant is that it measures exactly how large the Lagrange basis can amplify nodal data. Let $g \in C([-1,1])$ be arbitrary. From the Lagrange representation,
\begin{align*}
(I_N g)(x)=\sum_{j=0}^N g(x_j)\ell_j(x)
\end{align*}
for every $x \in [-1,1]$.
We now estimate this expression using only the supremum norm of $g$. For a fixed $x \in [-1,1]$, the triangle inequality gives
\begin{align*}
|(I_N g)(x)|\leq \sum_{j=0}^N |g(x_j)| |\ell_j(x)|.
\end{align*}
Each node $x_j$ lies in $[-1,1]$, so the definition of the supremum norm gives $|g(x_j)|\leq \|g\|_\infty$. Substituting this bound into the previous inequality yields
\begin{align*}
|(I_N g)(x)|\leq \|g\|_\infty \sum_{j=0}^N |\ell_j(x)|.
\end{align*}
This holds for every $x \in [-1,1]$. Therefore, taking the maximum over $x$ and using the definition
\begin{align*}
\lambda_N=\max_{x \in [-1,1]} \sum_{j=0}^N |\ell_j(x)|
\end{align*}
gives
\begin{align*}
\|I_N g\|_\infty \leq \lambda_N \|g\|_\infty.
\end{align*}
So the interpolation operator cannot increase the uniform norm by more than the Lebesgue constant.
[/guided]
[/step]
[step:Decompose the interpolation error using an arbitrary comparison polynomial]
Fix $p \in \mathcal{P}_N$. Since $I_N p=p$, linearity of $I_N$ gives
\begin{align*}
f-I_N f=f-p-I_N(f-p).
\end{align*}
Taking supremum norms and applying the triangle inequality gives
\begin{align*}
\|f-I_N f\|_\infty \leq \|f-p\|_\infty+\|I_N(f-p)\|_\infty.
\end{align*}
Using the bound from the previous step with $g=f-p$ gives
\begin{align*}
\|I_N(f-p)\|_\infty \leq \lambda_N\|f-p\|_\infty.
\end{align*}
Therefore
\begin{align*}
\|f-I_N f\|_\infty \leq (1+\lambda_N)\|f-p\|_\infty.
\end{align*}
[/step]
[step:Take the best polynomial approximation error]
The preceding estimate holds for every $p \in \mathcal{P}_N$. Hence taking the infimum over all such $p$ yields
\begin{align*}
\|f-I_N f\|_\infty \leq (1+\lambda_N)\inf_{p\in\mathcal P_N}\|f-p\|_\infty.
\end{align*}
This is the claimed Lebesgue inequality.
[/step]
