Recovering original regression coefficients from standardized Suppose I use Least Squares to estimate coefficients in the standard linear model with design matrix $X$'s columns standardized, so the model is
$$
E[y] = X^*\beta^*
$$
where $X^*$ is $X$ with columns centered and scaled. 
Assume $X$ has full column rank. To recover the regression coefficients $\beta$ in the model with unstandardized X
$$
E[y] = X\beta
$$
I should be able to use the equation
$$
\qquad \quad X\beta = X^*\beta^* \\
\implies X^T X\beta = X^TX^*\beta^* \\
\qquad \quad \; \; \; \; \implies \qquad \beta = (X^T X)^{-1}X^TX^*\beta^*
$$
To test this I ran some simple R code:  
set.seed(100)
# set number of samples
N <- 10
# set number of regressors
p <- 3
y <- runif(N)
X <- matrix(runif(N*p), N)
# Least squares coefficients with unstandardized X
(b1 <- solve(crossprod(X), crossprod(X,y)))

# now standardize X and get new coefficients
Xs <- scale(X)
(b2 <- solve(crossprod(Xs), crossprod(Xs,y)))

# b.orig should be exactly the same as b1, but its not!
(b.orig <- solve(crossprod(X), crossprod(X, Xs %*% b2)))

Here is the output of b1:  
            [,1]
[1,]  0.17109189
[2,]  0.52204169
[3,] -0.02115178

and the output of b.orig, which should equal b1, but does not:  
            [,1]
[1,] -0.15935376
[2,]  0.16915433
[3,] -0.04696604

What is going wrong here? I have already looked at:http://bit.ly/1HxHS6U, who seems to use a similar idea (equating expected or fitted values of y in both equations). 
 A: The models are not the same.  Therefore the coefficients should differ.  When you recenter, you are taking linear combinations of the columns of $X$ with the vector $\mathbf{1}=(1,1,\ldots, 1)^\prime$.  This is fine, provided that $\mathbf{1}$ lies in the column space.  In your example it does not.
What is worse, when you do include $\mathbf{1}$ as a column, the entire calculation falls apart due to singularities.
In detail, the original model (including a constant term) should be
$$\mathbb{E}(Y) = \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p = \beta_0 + \sum_{i=1}^p \beta_i X_i.$$
Standardizing to $X_i = \sigma_i Z_i + \mu_i$ yields
$$\mathbb{E}(Y) = \beta_0 + \sum_{i=1}^p \beta_i (\sigma_i Z_i + \mu_i) = \left(\beta_0 + \sum_{i=1}^p \beta_i \mu_i\right) + \sum_{i=1}^p (\beta_i \sigma_i) Z_i = \beta_0^{*} + \sum_{i=1}^p \beta_i^{*} Z_i,$$
with $\beta_i^{*} = \sigma_i \beta_i$ giving the correct relationships between the "standardized" and unstandardized coefficients.  (The usual definition of standardized coefficient also standardizes the response variable, so these $\beta_i^{*}$ might better be called "semi-standardized.")
