CD least squares

gistfile1.txt

import numpy as np
from numba import njit
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz
from numpy.linalg import norm
from sklearn.utils import check_random_state
from sklearn.linear_model import LinearRegression


def data(n_samples, n_features, rho=0.5, seed=24):
    corr = rho ** np.arange(n_features)
    cov = toeplitz(corr)

    rng = check_random_state(seed)
    X = rng.multivariate_normal(np.zeros(n_features), cov, n_samples)
    return np.asfortranarray(X)


@njit
def cd(X, y, max_iter):
    n_features = X.shape[1]
    R = y.copy()
    w = np.zeros(n_features)
    lc = (X ** 2).sum(axis=0)
    E = []
    for t in range(max_iter):
        for j in range(n_features):
            old_wj = w[j]
            w[j] += X[:, j].T @ R / lc[j]
            R += (old_wj - w[j]) * X[:, j]
        E.append((R @ R) / 2.)

    return w, np.array(E)


if __name__ == "__main__":
    n_samples = 100
    y = np.random.randn(n_samples)
    plt.close('all')
    fig, axarr = plt.subplots(2, 2, constrained_layout=True)

    for i, n_features in enumerate([50, 99, 100, 110]):
        X = data(n_samples, n_features, rho=0.9)

        w, E = cd(X, y, 200000)

        clf = LinearRegression(fit_intercept=False)
        clf.fit(X, y)

        E_star = norm(y - X @ clf.coef_) ** 2 / 2.

        good = np.where((E - E_star) > 1e-14)[0]
        axarr.flat[i].semilogy(good, (E - E_star)[good])
        axarr.flat[i].set_title("n_features = % d" % n_features)

    plt.show(block=False)