Confusion about ensembles and averages in autocorrelation matrices
EDIT: sorry that I didn't fully get your question. Now I got your points of conufion and try to provide a shortest possible answer:

First, you should better use xcorr function when estimating ACS $r_x[k]$ from data sequence $x[n]$. The output of xcorr will return the lag-0 at the middle sample as you recognize. Look at the input argument scaleopt which provides normalization to the computed output. The default is none ! and causes a confusion for you. You should select biased so that equation 7 below and xcorr output will be identical. Like this:

rx = xcorr([1,2,3],'biased'); % identical to (1/3)*xcorr([1,2,3])


Second,  do not use the raw output of xcorr as the input to the toeplitz function. As the output of xcorr will include lags for rx[-2],rx[-1],rx[0],rx[1],rx[2]. But you only need rx[0],rx[1],rx[2]. So you better call it like:

L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) );  % now yo uare ok ;-) 


Now it could be argued that the Matlab's default xcorr scale option none is misleading. And yes, I would say so. In its default form it computes the convolution and leaves it unscaled...

Finally, please distingusih a theoretical expectation from practical estimation. When algorithms work on data, they estimate and they don't compute expectations, which is a theoretical tool. Hence, all expectations are replaced with estimations when algorithms are run for the practical applications...
It seems (to me) that your confusion is about three things.


Theoretical definition of the autocorrelation sequence $r_{XX}[k]$ of a random process $X[n,s]$
Practical estimation of it from a sample function $x[n]$ via ergodicity.  
Interpretation of a random vector (including samples of a random process)


Consider the interpretation of a random process $X[n,s]$ as an ordered collection (set) of indexed random variables $X_n(s)$ for each time index $n$. Each random variable $X_n$ may, in principle, has its own unique probability distribution/density function $F_{X_n}(x)$. Consider the enormous implication of this on the number and type of all possible joint density functions.

To make things extremely simplified, we make the assumptions that all those random variables have the same PDF $F_{X_n}(x)$ for all $n$, and independent of each other, also further that first and second moments are independent of partciular time instants, which yields an i.i.d. , WSS random process $X[n,s]$ at hand.

The theoretical auto-correlation sequence (ACS) of this WSS random process is defined as:

$$ r_{X_n X_m}(n,m) = E\{X_n X^*_m\} = E\{X[n]X^*[n+k] \} = r_{X X}[k]   \tag{1}$$ 

The computation of this theoretical ACS $r_x[k]$ requires you to evaluate the integrals with PDF functions. This is possible with questions on the paper, but with applications in working systems, where you won't solve (mathematical) questions but perform statistical calculations, you need to estimate it from data, based on ergodicty assumption, which states that theoretical ensemble averages can bu estimates from practical sample function time-averages. Then you have

$$ \hat{r}_X[k] = \lim_{N \to \infty} \frac{1}{2N+1} \sum_{n=-N}^{N} x[n]x^*[n+k] \tag{2}$$

Note that $x[n]$ is a particular sample function (a realization) of the hrandom process $X[n,s]$. Note also that the infinite limits are replaced by data observation boundaries in practice.

Now consider the random vector $X(n) = [X(n),X(n-1),..,X(n-N+1)]^T$ which includes $N$ consecutive random variables from the random process $X[n,s]$, beginning at index $n$.

From this vector build up the outer product $X \cdot X^H$ , which includes a matrix of pairwise multiplications of random variables $X_i X_j$ as in

$$
RM =
  \begin{bmatrix}
X_1 X^*_1 &amp; X_1 X^*_2 &amp; ... &amp; X_1 X^*_N \\
X_2 X^*_1 &amp; X_2 X^*_2 &amp; ... &amp; X_2 X^*_N \\
... &amp; ... &amp; ... &amp; ... \\
X_N X^*_1 &amp; X_N X^*_2 &amp; ... &amp; X_N X^*_N \\
  \end{bmatrix} \tag{3}
$$

Now taking the expectation of this $RM$ matrix yields the theoretical auto-correlation matrix of the WSS random process as:

$$ R_X = E\{ X \cdot X^H \} = E\{RM\} = E\{\begin{bmatrix}
X_1 X^*_1 &amp; X_1 X^*_2 &amp; ... &amp; X_1 X^*_N \\
X_2 X^*_1 &amp; X_2 X^*_2 &amp; ... &amp; X_2 X^*_N \\
... &amp; ... &amp; ... &amp; ... \\
X_N X^*_1 &amp; X_N X^*_2 &amp; ... &amp; X_N X^*_N \\
  \end{bmatrix} \} \tag{4}$$

and this results in:
$$ R_X = E\{ X \cdot X^H \} = \begin{bmatrix}
r_x[0] &amp; r_x[-1] &amp; ... &amp; r_x[-N+1] \\
r_x[1] &amp; r_x[0] &amp; ... &amp; r_x[-N+2] \\
... &amp; ... &amp; ... &amp; ... \\
r_x[N-1] &amp; r_x[N-2] &amp; ... &amp; r_x[0] \\
  \end{bmatrix} \tag{5}
 $$

for real data it becomes:
$$ R _x = 
 \begin{bmatrix}
r_x[0] &amp; r_x[1]^* &amp; ... &amp; r_x[N-1]^* \\
r_x[1] &amp; r_x[0] &amp; ... &amp; r_x[N-2]^* \\
... &amp; ... &amp; ... &amp; ... \\
r_x[N-1] &amp; r_x[N-2] &amp; ... &amp; r_x[0] \\
  \end{bmatrix} \tag{6}
$$

Again this theoretical ACM $R_X$ cannot be computed but is estimated from sample function data...

Estimation of the auto-correlation matrix $R_x$ can be done in a number of ways (both related) but at the basic level you see that what you should actually estimate is the auto-correlations sequence $r_x[k]$ for the lags $k=0,1,...,N-1$. Then from the Toeplitz and Hermitian symmetric property of $R_x$, you can construct it. 

One way to practically estimate ACS $r_x[k]$ from the data $x[n]$ of length $L$ is :

$$ \hat{r}[k] = \frac{1}{L} \left( x[k] \star x^*[-k] \right) \tag{7}$$ 

This yields an estimate of ACS for the lags $-L &lt; k  &lt; L$, (you only need $N$ of those for $k=0,1,...,N-1$)

Alternatively, ACM can be estimated as a whole similarly to ACS from: 

$$ \hat{R}_x = \frac{1}{K} \sum_{k=1}^{K} X(k) \cdot X(k)^H \tag{8}$$  where data vector $X(k)$ includes $N$ consecutive samples from $x[n]$ (specially processed at the edges)