In the proof of Theorem 2 of Sparse projections onto the simplex authors mention the following equality for any $\mathbf{b} \in \mathbb{R}^k$ and $\lambda \in \mathbb{R}$:
$$ \begin{aligned} &\sum_{i=1}^k b_i^2 - \frac{((\sum_{i=1}^k b_i) - \lambda)^2}{k} \\ &= \\ &\lambda(2b_1-\lambda)+\sum_{j=2}^k \frac{j-1}{j} \Big( b_j - \frac{(\sum_{i=1}^{j-1}b_i)-\lambda}{j-1} \Big)^2. \end{aligned} $$
Question
Is there an elegant way of showing the above equality? Probably it is possible to show it using induction but I am looking to find a way other than induction.
My thoughts
If we look at the way terms are given on the left hand side, one can think of sampling $k$ number of $b_i$'s with uniform probability of $\frac{1}{k}$ where each $b_i$ has a fix distribution $\mathcal{D}$. Then the sums can be written as expected value of k samples of $b_i$ and $b_i^2$. Looking in this way, can we give another proof to that other than algebraic proof?
Another observation is that, the right hand side some how tries to relate sum of squares and square of sums in a recursive way.
Seeing LHS as a quadratic expression of $\lambda$, LHS can be written as $$A\lambda^2+B\lambda+C$$ where $$A=-\dfrac 1k,\qquad B=\dfrac{2}{k}\sum_{i=1}^{k}b_i,\qquad C=\bigg(\sum_{i=1}^{k}b_i^2\bigg)-\dfrac 1k\bigg(\sum_{i=1}^{k}b_i\bigg)^2$$
Seeing RHS as a quadratic expression of $\lambda$, RHS can be written as $$D\lambda^2+E\lambda+F$$ where $$D=-1+\sum_{j=2}^{k}\frac{1}{j(j-1)}$$ $$E=2b_1+\bigg(\sum_{j=2}^{k}\frac{2b_j}{j}\bigg)-\sum_{j=2}^{k}\frac{2}{j(j-1)}\sum_{i=1}^{j-1}b_i$$ $$F=\bigg(\sum_{j=2}^{k}\frac{j-1}{j}b_j^2\bigg)-\bigg(\sum_{j=2}^{k}\frac{2b_j}{j}\sum_{i=1}^{j-1}b_i\bigg)+\sum_{j=2}^{k}\frac{1}{j(j-1)}\bigg(\sum_{i=1}^{j-1}b_i\bigg)^2$$
In the following, let us prove that $A=D, B=E$ and $C=F$ using the idea of telescoping sum : $$\sum_{j=a}^{b}\frac{1}{j(j-1)}=\sum_{j=a}^{b}\bigg(\frac{1}{j-1}-\frac 1j\bigg)=\frac{1}{a-1}-\frac 1b$$
Proof for $A=D$ :
We have $$D=-1+\sum_{j=2}^{k}\frac{1}{j(j-1)}=-1+\bigg(1-\frac 1k\bigg)=-\frac 1k=A$$
Proof for $B=E$ :
$(\text{coefficient of $b_1$ in $E$})=2-2\displaystyle\sum_{j=2}^{k}\dfrac{1}{j(j-1)}=\dfrac 2k$
$(\text{coefficient of $b_m$ in $E$ where $1\lt m\lt k$})=\dfrac{2}{m}-2\displaystyle\sum_{j=m+1}^{k}\dfrac{1}{j(j-1)}=\dfrac 2k$
$(\text{coefficient of $b_k$ in $E$)}=\dfrac 2k$
It follows from these that $B=E$.
Proof for $C=F$ :
$(\text{coefficient of $b_1^2$ in $F$})=\displaystyle\sum_{j=2}^{k}\dfrac{1}{j(j-1)}=1-\dfrac 1k$
$(\text{coefficient of $b_m^2$ in $F$ where $1\lt m\lt k$})=\dfrac{m-1}{m}+\displaystyle\sum_{j=m+1}^{k}\dfrac{1}{j(j-1)}=1-\dfrac 1k$
$(\text{coefficient of $b_k^2$ in $F$)}=1-\dfrac 1k$
$(\text{coefficient of $b_mb_n$ in $F$ where $m>n$})=-\dfrac{2}{m}+2\displaystyle\sum_{j=m+1}^{k}\dfrac{1}{j(j-1)}=-\dfrac 2k$
It follows from these that $C=F$.
Since $A=D,B=E$ and $C=F$, we see that $\text{LHS $=$ RHS}$ does hold.$\ \blacksquare$
Added :
We have $$\begin{align}F&=\bigg(\sum_{j=2}^{k}\frac{j-1}{j}b_j^2\bigg)-\bigg(\sum_{j=2}^{k}\frac{2b_j}{j}\sum_{i=1}^{j-1}b_i\bigg)+\sum_{j=2}^{k}\frac{1}{j(j-1)}\bigg(\sum_{i=1}^{j-1}b_i\bigg)^2 \\\\&=\frac{2-1}{2}b_2^2+\frac{3-1}{3}b_3^2+\cdots +\frac{k-1}{k}b_k^2 \\\\&\quad -\frac{2b_2}{2}b_1-\frac{2b_3}{3}(b_1+b_2)-\cdots -\frac{2b_k}{k}(b_1+b_2+\cdots +b_{k-1}) \\\\&\quad +\frac{b_1^2}{2(2-1)}+\frac{(b_1+b_2)^2}{3(3-1)}+\cdots +\frac{(b_1+b_2+\cdots +b_{k-1})^2}{k(k-1)}\end{align}$$