Is this a mistake in the feature standardization of ridge regression in book Machine Learning in Action?

Question

I am new to ML. I am currently reading the classic book Machine Learning in Action by Peter Harrington. In its implementation of ridge regression in P165 Listing 8.3, the book standardizes feature matrix $X$ by subtracting means of attributes and then dividing it by variances of attributes not by standard deviations of attributes! As follows:

def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat
    denom = xTx + eye(shape(xMat)[1])*lam
    if linalg.det(denom) == 0.0:
        print "This matrix is singular, cannot do inverse"
        return
    ws = denom.I * (xMat.T*yMat)
    return ws

def ridgeTest(xArr,yArr):
    xMat = mat(xArr); yMat=mat(yArr).T
    yMean = mean(yMat,0)
    yMat = yMat - yMean     #to eliminate X0 take mean off of Y
    #regularize X's
    xMeans = mean(xMat,0)   #calc mean then subtract it off
    xVar = var(xMat,0)      #calc variance of Xi then divide by it
    xMat = (xMat - xMeans)/xVar
    numTestPts = 30
    wMat = zeros((numTestPts,shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,exp(i-10))
        wMat[i,:]=ws.T
    return wMat

Then calls ridgeTest by

>>> abX,abY=regression.loadDataSet('abalone.txt')
>>> ridgeWeights=regression.ridgeTest(abX,abY)

I do not think this makes any sense. According to standardization, To create a unit variance, we should divide $X-\mu$ by standard deviation $\sigma$, because: $$ D(\frac{X-\mu}{\sigma})=\frac{1}{\sigma^2}D(X)=\frac{1}{\sigma^2}\sigma^2=1, $$ where $D(X)=\sigma^2$ is the variance of $X$.

If according to the book, divide $X-\mu$ by variance $\sigma^2$, we only get: $$ D(\frac{X-\mu}{\sigma^2})=\frac{1}{\sigma^4}D(X)=\frac{1}{\sigma^4}\sigma^2=\frac{1}{\sigma^2}, $$ which is not a unit variance.

The reason why I am not so confident that this is a mistake is that, first of all, the book is handling $Y$ carefully by subtracting its mean and not include an all 1 column in $X$ passed, both of which I find to be proper (see L2-normalization does not punish intercept). Hence I do not believe the book will make an obvious mistake like this. Secondly, nor the errata (see errata) or the book forum discussed it before. It is not a new book.

Although I posted divide by std not variance in the book forum, no reply has been received yet. So I turn to stack exchange. Does the feature standardization of the book mistakenly divide $X-\mu$ by variance or it is a meaningful action?

Thanks a lot in advance! Any clue will help me.

Answer 1

Short answer

Yes, it seems to be a mistake and may have an imporant impact on the way some variables are penalized.

If you have two variables with different magnitude, this rescaling with basically give a low standard deviation to the value with the (initial) largest standard deviation. Therefore, this variable will suffer a smaller penalization than it should. And this will make your regression depend on the unit in which the values are expressed (which is not a desirable feature).

Side note

The algorithm scales the input, but returns the outputs evaluated on the scaled inputs. This makes it quite hard to implement a "predict" function...

Numerical experiment

I fixed the issue (in the case where there is not intercept) so that the algorithm returns wMat[i,:]=ws.T / xVar.

import numpy as np
import matplotlib.pyplot as plt

def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat
    denom = xTx + np.eye(np.shape(xMat)[1])*lam
    if np.linalg.det(denom) == 0.0:
        print "This matrix is singular, cannot do inverse"
        return
    ws = denom.I * (xMat.T*yMat)
    return ws

# Assuming there is not intercept : y = x_0 + 2 * x_1

def ridgeTest(xArr,yArr):
    xMat = np.mat(xArr); yMat=np.mat(yArr).T
    yMean = np.mean(yMat,0)
    yMat = yMat - yMean     #to eliminate X0 take mean off of Y
    #regularize X's
    xMeans = np.mean(xMat,0)   #calc mean then subtract it off
    xVar = np.var(xMat,0)      #calc variance of Xi then divide by it
    xMat = (xMat - xMeans)/xVar
    numTestPts = 30
    wMat = np.zeros((numTestPts,np.shape(xMat)[1]))
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,np.exp(i-10))
        wMat[i,:]=ws.T / xVar
    return wMat 

X = np.mat(np.random.rand(50,2))
y = np.mat([1,2]) * X.T 

res_orig = ridgeTest(X,y)

plt.plot(res_orig)
plt.title('Original Scaling')
plt.ylabel('Values of the (estimated) coefficient.')
plt.xlabel('Value of the penalization.')
plt.show()


X[:,1] *= 10
y = np.mat([1,2]) * X.T 
res = ridgeTest(X,y)

plt.plot(res)
plt.title('x_1 has been multiplied by 10')
plt.ylabel('Values of the (estimated) coefficient.')
plt.xlabel('Value of the penalization.')
plt.show()

Result

The result with the original scaling provide the actual values ($y=x_0+2x_1$) whereas the second case exhibits a strange behavior.

Edit

As stated by @Berkmeister in the comments, this scaling can be defended : in many models, you can decide, arbitrarly, to change the weights of your features.

But, the fact that the coefficients returned by the model are not rescaled, and that this choice of (uncommon) scaling is proposed without prior explanations let me think that this code can be improved.