Shrinkage of the Sample Covariance matrix Assume we have N independent and identically distributed random vectors $X_1, X_2, ..., X_N$ where each of them is of size p $\times$ 1. The sample covariance matrix, denoted here by $S$, is computed as $(1/N)XX'$ where $X$ is now a matrix having the $X_i$ as columns, $i = 1, 2, ..., N$.
As we all know that $S$, has no structure, unbiased but has a lot of estimation error(ill-conditioned). In addition, $S$ becomes singular when $p>=N$. 
So the obvious idea is to try by shrinking the unbiased covariance matrix $S$ towards a biased (with less variance) structured covariance often called "The shrinkage Target". The shrinkage formula can be expressed as follows:
$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$$(1-\alpha)S + \alpha F$ 
Where $\alpha$ is the shrinkage intensity (between 0 and 1), and $F$ is the shrinkage target.
In shrinkage, the most challenging part is how to compute (automatically) the shrinkage intensity. So following the article of Ledoit and Wolf in 2003 "Improved Estimation of the Covariance Matrix of Stock Returns
With an Application to Portfolio Selection", they computed an automatic optimal shrinkage intensity after minimizing a certain loss function based on the Frobenius distance.
Ok great! My question now is from this article. In page 8, I want to know how they just computed $[E(\alpha f_{ij} + (1 - \alpha)s_{ij} - \delta_{ij}]^2$. Ok, in the article, they used the single index model as shrinkage target $F$. But in fact, let us compute this squared expected value for a general shrinkage target. I just want to know how they computed this squared expectation (I am just interested for this expectation in page 8 before they did the derivation w.r.t $\alpha$)...
Kindly any help will be very very appreciated!! 
 A: What they evaluate there is
$$ \mathbb{E}\big[ (\alpha f + (1-\alpha)s - \sigma)^2 \big], $$
(note this isn't identical to what you wrote in the question) where you are given
$$ \mathbb{E}[f] = \phi, \qquad \mathbb{E}[s] = \sigma. $$
They evaluate this by writing
$$ \alpha f + (1-\alpha)s-\sigma =: X, $$
and using the identity
$$ \mathbb{E}[X^2] = \mathbb{V}[X] + (\mathbb{E}[X])^2. $$
Since $\phi$ and $\sigma$ are constants here,
$$ \mathbb{V}[X] = \alpha^2 \mathbb{V}[f] + (1-\alpha)^2 \mathbb{V}[s] + \alpha(1-\alpha)\mathrm{Cov}[f,s], $$
and 
$$ \mathbb{E}[X] = \alpha(\phi-\sigma) $$
due to $\mathbb{E}[f]=\phi$, $\mathbb{E}[s]=\sigma$, and the usual variance-of-sum formula.
To generalize this to any other $\phi$, note that none of the properties of their matrices $\mathbf{F}$, $\mathbf{S}$, $\mathbf{\Sigma}$ are actually being used: this is simply an identity in terms of random variables $f,s$ and their means $\phi,\sigma$. It does not matter whether $\phi$ are the identity matrix elements or something else.
