Regularised regression code

Question

I wrote some code to do regularised linear regression and it works, but I don't like the fact that I've had to double call the functions when plotting, nor the fact that I've sliced those calls to get the parts that I want.

import numpy as np
import matplotlib.pyplot as plt

def fit(phi_fn, xx, yy):
    w_fit = np.linalg.lstsq(phi_fn(xx), yy, rcond=None)[0]
    grid_size = 0.01
    x_grid = np.arange(0,9,grid_size)[:,None]
    f_grid = np.matmul(phi_fn(x_grid),w_fit)
    return(x_grid, f_grid)

def fitreg(phi_fn, xx, yy,lamb):
    yy = np.pad(yy,(0,8),'constant',constant_values=0)
    zz = np.concatenate((phi_fn(xx),lamb*np.identity(8)),axis=0)
    w_fit = np.linalg.lstsq(zz, yy, rcond=None)[0]
    grid_size = 0.01
    x_grid = np.arange(0,9,grid_size)[:,None]
    f_grid = np.matmul(phi_fn(x_grid),w_fit)
    return(x_grid, f_grid)

def phi_poly(xx):
        return np.concatenate([np.ones((xx.shape[0],1)), xx,xx**2,xx**3,xx**4,xx**5,xx**6,xx**7], axis=1)


D = 1
N = 10
mu = np.array([0,1,2,3,4,5,6,7,8,9])

xx = np.tile(mu[:,None], (1, D)) + 0.01*np.random.randn(N, D)
yy = 2*xx + 2*np.random.randn(N,D)

plt.clf()
plt.plot(xx,yy,'kx')

plt.plot(fit(phi_poly, xx, yy)[0], fit(phi_poly, xx, yy)[1], 'b-')
plt.plot(fitreg(phi_poly, xx, yy,1)[0], fitreg(phi_poly, xx, yy,1)[1][:,0], 'r-')
plt.plot(fitreg(phi_poly, xx, yy,10)[0], fitreg(phi_poly, xx, yy,10)[1][:,0], 'g-')
plt.plot(fitreg(phi_poly, xx, yy,0.1)[0], fitreg(phi_poly, xx, yy,0.1)[1][:,0], 'y-')

plt.show()

Answer 1

You can get rid of some of the repetition using a for loop looping over the lambdas and colors and using tuple assignment to at least get rid of one level of having to use indexing as well as the necessity of having to call fitreg twice for each plot:

if __name__ == "__main__":

    D = 1
    N = 10
    mu = np.arange(10)

    xx = mu[:, None] + 0.01 * np.random.randn(N, D)
    yy = 2 * xx + 2 * np.random.randn(N, D)

    plt.clf()
    plt.plot(xx, yy, 'kx')

    x, y = fit(phi_poly, xx, yy)
    plt.plot(x, y, 'b-')
    for lamb, color in zip([1, 10, 0.1], "rgy"):
        x, y = fitreg(phi_poly, xx, yy, lamb)
        plt.plot(x, y[:, 0], '{}-'.format(color))

    plt.show()

I also added a if __name__ == "__main__": guard to allow importing from this script in another script and used numpy.arange for your mu array.

For the construction of xx you don't actually need the numpy.tile, the normal broadcasting rules already do that for you.