Where is the explanatory effect of common variance among covariates accounted for in regression procedures? As a follow up to the excellent answers provided for:
Does the order of explanatory variables matter when calculating their regression coefficients?
(Which I've found incredibly useful from a pedagogical perspective), I've been wondering how exactly it manages to provide regression coefficients when we deal with high collinear data (setting aside the high standard error of these 
estimates).
Edit: For ease I've reproduced the section in the linked question that gets to the crux of the confusion (from Elements of Statistical Learning), the first two image provide the background but the section in italics in the final image gets to the root of the intuition which I am struggling with:



My question in words, is - if, as stated above, multiple regression coefficients express the effect of each covariate on a dependent variable having partialled out the variability that can be explained by other variables, where is the explanatory effect of the shared variability of the covariates accounted for?
Note I am hoping to get the intuition here - the algebra and geometry of the solution and both fairly easy to grasp.
As an example that tries to elucidate, consider a logical extreme where:
$$ Y = X + \epsilon_y $$
$$ \epsilon_y \sim N(0,0.1) $$ 
$$ X_1 = X + \epsilon_1 $$
$$ X_2 = X + \epsilon_2 $$
$$ \epsilon_1 \sim \epsilon_2 \sim N(0,0.001) $$ 
That is, $Y$ and $X$ have a strong linear relationship, and there is strong collinearity between $X_1$ and $X_2$ caused by their common factor $X$. Now suppose we attempt:
$$ Y \sim X_1 + X_2 $$
Following the Gram-Schmidt procedure, the residual of $X_1$ or $X_2$ on the other covariates (in this case, just each other) effectively remove the common variance between them (this may be where I am misunderstanding), but surely doing so removes the common element that manages to explain the relationship with $Y$?
Edit: To clarify on a point made by below: as is elaborated in the linked question, in the GS procedure the multiple regression coefficients are not generated from the interim coefficients that are produced 'en route' to the final residual. That is, to get the coefficient for $X_2$ we take the GS procedure from intercept > $X_1$ > $X_2$. Then to generate the coefficient for $X_1$ we would work through intercept > $X_2$ > $X_1$. In both instances the crucial common variance due to $X$ and the resulting relationship with $Y$ is lost. 
 A: Even though you say that the geometry of this is fairly clear to you, I think it is a good idea to review it. I made this back of an envelope sketch:

Left subplot is the same figure as in the book: consider two predictors $x_1$ and $x_2$; as vectors, $\mathbf x_1$ and $\mathbf x_2$ span a plane in the $n$-dimensional space, and $\mathbf y$ is being projected onto this plane resulting in the $\hat {\mathbf y}$.
Middle subplot shows the $X$ plane in the case when $\mathbf x_1$ and $\mathbf x_2$ are not orthogonal, but both have unit length. The regression coefficients $\beta_1$ and $\beta_2$ can be obtained by a non-orthogonal projection of $\hat{\mathbf  y}$ onto $\mathbf x_1$ and $\mathbf x_2$: that should be pretty clear from the picture. But what happens when we follow the orthogonalization route?
The two orthogonalized vectors $\mathbf z_1$ and $\mathbf z_2$ from Algorithm 3.1 are also shown on the figure. Note that each of them is obtained via a separate Gram-Schmidt orthogonalization procedure (separate run of Algorithm 3.1): $\mathbf z_1$ is the residual of $\mathbf x_1$ when regressed on $\mathbf x_2$ ans $\mathbf z_2$ is the residual of $\mathbf x_2$ when regressed on $\mathbf x_1$. Therefore $\mathbf z_1$ and $\mathbf z_2$ are orthogonal to $\mathbf x_2$ and $\mathbf x_1$ respectively, and their lengths are less than $1$. This is crucial.
As stated in the book, the regression coefficient $\beta_i$ can be obtained as $$\beta_i = \frac{\mathbf z_i \cdot \mathbf y}{\|\mathbf z_i\|^2} =\frac{\mathbf e_{\mathbf z_i} \cdot \mathbf y}{\|\mathbf z_i\|},$$ where $\mathbf e_{\mathbf z_{i}}$ denotes a unit vector in the direction of $\mathbf z_i$. When I project $\hat{\mathbf  y}$ onto $\mathbf z_i$ on my drawing, the length of the projection (shown on the figure) is the nominator of this fraction. To get the actual $\beta_i$ value, one needs to divide by the length of $\mathbf z_i$ which is smaller than $1$, i.e. the $\beta_i$ will be larger than the length of the projection.
Now consider what happens in the extreme case of very high correlation (right subplot). Both $\beta_i$ are sizeable, but both  $\mathbf z_i$ vectors are tiny, and the projections of $\hat{\mathbf  y}$ onto the directions of $\mathbf z_i$ will also be tiny;  this is I think what is ultimately worrying you. However, to get $\beta_i$ values, we will have to rescale these projections by inverse lengths of $\mathbf z_i$, obtaining the correct values.

Following the Gram-Schmidt procedure, the residual of X1 or X2 on the other covariates (in this case, just each other) effectively remove the common variance between them (this may be where I am misunderstanding), but surely doing so removes the common element that manages to explain the relationship with Y?

To repeat: yes, the "common variance" is almost (but not entirely) "removed" from the residuals -- that's why projections on $\mathbf z_1$ and $\mathbf z_2$ will be so short. However, the Gram-Schmidt procedure can account for it by normalizing by the lengths of $\mathbf z_1$ and $\mathbf z_2$; the lengths are inversely related to the correlation between $\mathbf x_1$ and $\mathbf x_2$, so in the end the balance gets restored.

Update 1
Following the discussion with @mpiktas in the comments: the above description is not how Gram-Schmidt procedure would usually be applied to compute regression coefficients. Instead of running Algorithm 3.1 many times (each time rearranging the sequence of predictors), one can obtain all regression coefficients from the single run. This is noted in Hastie et al. on the next page (page 55) and is the content of Exercise 3.4. But as I understood OP's question, it referred to the multiple-runs approach (that yields explicit formulas for $\beta_i$).
Update 2
In reply to OP's comment:

I am trying to understand how 'common explanatory power' of a (sub)set of covariates is 'spread between' the coefficient estimates of those covariates. I think the explanation lies somewhere between the geometric illustration you have provided and mpiktas point about how the coefficients should sum to the regression coefficient of the common factor

I think if you are trying to understand how the "shared part" of the predictors is being represented in the regression coefficients, then you do not need to think about Gram-Schmidt at all. Yes, it will be "spread out" between the predictors. Perhaps a more useful way to think about it is in terms of transforming the predictors with PCA to get orthogonal predictors. In your example there will be a large first principal component with almost equal weights for $x_1$ and $x_2$. So the corresponding regression coefficient will have to be "split" between $x_1$ and $x_2$ in equal proportions. The second principal component will be small and $\mathbf y$ will be almost orthogonal to it.
In my answer above I assumed that you are specifically confused about Gram-Schmidt procedure and the resulting formula for $\beta_i$ in terms of $z_i$.
A: The GS procedure would start with $X_1$ and then move to orthogonalizing $X_2$. Since $X_1$ and $X_2$ share $X$ the result would practically be zero in your example. But the common element $X$ remains, because we started with $X_1$, and $X_1$ still has $X$. 
Since $X_1$ and $X_2$ share common $X$, we would get that the remainder of $X_2$ after orthogonalization is practically zero as is stated in citation. 
In this case one could argue that original multiple regression problem is ill posed, so there is no sense to proceed, i.e. we should stop GS process and restate original multiple regression problem as $Y\sim X_1$. In this case we do not lose the common factor $X$ and correctly disregard $X_2$, since it does not give us any new information which we do not have. 
Of course we can proceed with GS procedure and calculate the coefficient for $X_2$ and recalculate back to the original multiple regression problem. Since we do not have perfect colinearity it is possible to do that theoretically. Practically it will depend on the numerical stability of the algorithms. Since 
$$\alpha X_1+ \beta X2 = (\alpha+\beta)X +\alpha\epsilon_1 + \beta\epsilon_2 $$
the regression  $Y\sim X_1 + X_2$ will produce coefficients $\alpha$ and $\beta$ such that $\alpha+\beta \approx 1$ (we will not have strict equality because of $\epsilon_1$ and $\epsilon_2$). 
Here is the example in R:
> set.seed(1001)
> x<-rnorm(1000)
> y<-x+rnorm(1000, sd = 0.1)
> x1 <- x + rnorm(1000, sd =0.001)
> x2 <- x + rnorm(1000, sd =0.001)
> lm(y~x1+x2)

Call:
lm(formula = y ~ x1 + x2)

Coefficients:
(Intercept)           x1           x2  
 -0.0003867   -1.9282079    2.9185409  

Here I skipped GS procedure, because the lm gave feasible results, and in that case recalculating coefficients from GS procedure does not fail.
