<jupyter_start><jupyter_text>Finite case$x\in[M]$<jupyter_code>M = 1000<jupyter_output><empty_output><jupyter_text>True function<jupyter_code>_d = 10
param = np.random.uniform(size=_d) - 1/2
fmap = np.array([[(x/M)**y for y in range(_d)] for x in range(M)])
#fmap = np.random.uniform(size = (M, _d))
_f = np.matmul(fmap, param)
f = lambda x: _f[x]
plt.plot(range(M), _f)
ylim = [min(_f)-0.1, max(_f)+0.1]<jupyter_output><empty_output><jupyter_text>Approximation with d random features (no noise)<jupyter_code>def rfmap(d):
    return np.random.normal(size=(M,d))
def ols(X, y, reg=0):
    A = np.matmul(X.T,X) + reg * np.eye(d)
    b = np.matmul(X.T, y)
    theta = np.linalg.solve(A,b)
    return np.matmul(X, theta), theta
RF = rfmap(M)
err = np.zeros(M)
perr = np.zeros(M)
fhats = np.zeros((M,M))
for d in range(1, M+1):
    X = RF[:, :d] 
    y = _f
    fhat, theta = ols(X, y)
    fhats[d-1] = fhat
    err[d-1] = np.linalg.norm(fhat - _f, ord=1) / M
    
plt.plot(range(M), err)
plt.plot(range(M), 0.2/np.sqrt(range(1,M+1)))
xx = range(M)

def update(d = 1):
    plt.plot(xx, _f, xx, fhats[d-1])
    axes = plt.gca()
    axes.set_ylim(ylim)
    
slider = widgets.interact(update, d=widgets.IntSlider(min=1, max=M, step=1, value=1))<jupyter_output><empty_output><jupyter_text>With noise<jupyter_code>N = 900
noise = 0.001
reg = 1
S = np.random.choice(M, size=N, replace=True)
y = _f[S] + noise * np.random.normal(size=N)
RF = rfmap(M)
err = np.zeros(M)
fhats = np.zeros((M,M))
for d in range(1, M+1):
    X = RF[S, :d] 
    fhat, theta = ols(X, y, reg)
    err[d-1] = np.linalg.norm(np.matmul(RF[:, :d], theta) - _f, ord=1) / N
    fhats[d-1] = np.matmul(RF[:, :d], theta)
plt.plot(range(M), err)
xx = range(M)

def update(d = 1):
    plt.scatter(S,y)
    plt.plot(xx, fhats[d-1])
    axes = plt.gca()
    axes.set_ylim(ylim)
    
slider = widgets.interact(update, d=widgets.IntSlider(min=1, max=M, step=1, value=M))<jupyter_output><empty_output><jupyter_text>Linear case$x\in\mathbb{R}^m$<jupyter_code>def make_fun(m):
    param = 2 * np.random.uniform(size=m) - 1
    return lambda x : np.matmul(x, param)
#This represents functions of the form <x,c> + d
def make_rfm(d, m):
    A = np.random.normal(size=(d, m))
    b = np.random.normal(size=d)
    return lambda X: np.matmul(A, X.T).T + b
def ls(A, b, reg=0):
    return np.linalg.solve(np.matmul(A.T, A) + reg * np.eye(A.shape[-1]), np.matmul(A.T, b))
def rf_approx(X, y, d, reg=0):
    rfm = make_rfm(d, X.shape[-1])
    Phi = rfm(X) #Nxd
    theta = ls(Phi, y, reg)
    yhat = np.matmul(Phi, theta)
    return yhat, rfm, theta
def gen_rank(n, m, r):
    A = 2 * np.random.uniform(size=(n, r)) - 1
    B = 2 * np.random.uniform(size=(r, r)) - 1
    C = 2 * np.random.uniform(size=(r, m)) - 1
    M = np.matmul(np.matmul(A, B), C)
    assert np.linalg.matrix_rank(M) == r
    return M<jupyter_output><empty_output><jupyter_text>Scalar example<jupyter_code>m = 1
N = 100
noise = 0
d = 2

f = make_fun(m)
X = 2 * np.random.uniform(size=(N, m)) - 1
y = f(X) + noise * np.random.normal(size=N)
yhat, _, _ = rf_approx(X, y, d)
xx = np.linspace(-1,1,100)
plt.scatter(X, y)
plt.plot(X, yhat)
ax = plt.gca()
ax.set_ylim([-1,1])
ax.set_xlim([-1,1])<jupyter_output><empty_output><jupyter_text>Full-rank case<jupyter_code>m = 100
N = 500
noise = 0.1
reg = 0
norm_ord = 1
max_d = 2 * m

f = make_fun(m)
X = 2 * np.random.uniform(size=(N, m)) - 1
y = f(X) + noise * np.random.normal(size=N)
ls_param = ls(X, y, reg)
yls = np.matmul(X, ls_param)
lserr = np.linalg.norm(yls - y, ord=norm_ord)**norm_ord / N
lsrank = np.linalg.matrix_rank(X)

X_test = 2 * np.random.uniform(size=(N, m)) - 1
y_test = f(X_test) + noise * np.random.normal(size=N)
yls_test = np.matmul(X_test, ls_param)
test_lserr = np.linalg.norm(yls_test - y_test, ord=norm_ord)**norm_ord / N

err = []
test_err = []
ranks = []
for d in range(1, max_d+1):
    yhat, rfm, theta = rf_approx(X, y, d, reg)
    err.append(np.linalg.norm(y - yhat, ord=norm_ord)**norm_ord / N)
    ranks.append(np.linalg.matrix_rank(rfm(X)))
    
    #Test
    yhat_test = np.matmul(rfm(X_test), theta)
    test_err.append(np.linalg.norm(y_test - yhat_test, ord=norm_ord)**norm_ord / N)

xx = range(1, max_d+1)
plt.plot(xx, err)
plt.plot(xx, test_err)
plt.plot([xx[0], xx[-1]], [lserr, lserr])
plt.plot([xx[0], xx[-1]], [test_lserr, test_lserr])
plt.xlabel('d')
plt.legend(('Training error', 'Test error', 'LS', 'LS test'))
xx = range(1, max_d+1)
plt.plot(xx, err)
plt.plot(xx, test_err)
plt.plot([xx[0], xx[-1]], [lserr, lserr])
plt.plot([xx[0], xx[-1]], [test_lserr, test_lserr])
plt.xlabel('d')
plt.legend(('Training error', 'Test error', 'LS', 'LS test'))
plt.gca().set_xlim((m-3, m+3))
plt.gca().set_ylim((0,err[m-4]))
print(lsrank, np.max(ranks))
plt.plot(xx, ranks)
plt.plot((xx[0], xx[-1]),(lsrank, lsrank))
plt.legend(('rank after rf', 'original rank'))
for i in range(max_d): print(i+1, err[i], test_err[i], '\n')<jupyter_output>1 2.488181510678169 2.4294863191972143 

2 2.4414629899031017 2.3851319361106413 

3 2.4670539805310914 2.4397782697215624 

4 2.3989992734110874 2.3523947999584607 

5 2.406287196946102 2.4491284244295546 

6 2.4258816131915886 2.3497609997047384 

7 2.4352556384967072 2.4258213169966356 

8 2.3980929116792913 2.3743782440824996 

9 2.386465232443141 2.3761536506451346 

10 2.316950626384221 2.314760042727451 

11 2.260506042904185 2.3532892699486756 

12 2.302896535642285 2.3700780665601116 

13 2.3086236456723173 2.3433184933934346 

14 2.2851281433587958 2.2925924062809804 

15 2.2419307862482873 2.3296093846056576 

16 2.3362047377801365 2.2798773634227087 

17 2.2038989626652947 2.2843289285009503 

18 2.092371606577256 2.1067207675538255 

19 2.260377188534007 2.2698325132127524 

20 2.1170364473096797 2.1642848092525426 

21 2.1370925601218604 2.1387564982718863 

22 2.1745914803313813 2.148056740130699 

23 2.146482777307317 2.2397467538113824 

24 2.2461750124274693 2.2426138[...]<jupyter_text>The rank-deficient case<jupyter_code>m = 100
N = 1000
noise = 0.1
reg = 0.1
norm_ord = 1
max_d = 3*m
rank = m - 1

f = make_fun(m)

#Contexts do not span R^d
X = gen_rank(N, m, rank)
lsrank = np.linalg.matrix_rank(X)

y = f(X) + noise * np.random.normal(size=N)
ls_param = ls(X, y, reg)
yls = np.matmul(X, ls_param)
lserr = np.linalg.norm(yls - y, ord=norm_ord)**norm_ord / N

X_test = gen_rank(N, m, rank)
y_test = f(X_test) + noise * np.random.normal(size=N)
yls_test = np.matmul(X_test, ls_param)
test_lserr = np.linalg.norm(yls_test - y_test, ord=norm_ord)**norm_ord / N

err = []
test_err = []
ranks = []
for d in range(1, max_d+1):
    yhat, rfm, theta = rf_approx(X, y, d, reg)
    err.append(np.linalg.norm(y - yhat, ord=norm_ord)**norm_ord / N)
    #print(rfm(X).shape)
    ranks.append(np.linalg.matrix_rank(rfm(X)))
    
    #Test
    yhat_test = np.matmul(rfm(X_test), theta)
    test_err.append(np.linalg.norm(y_test - yhat_test, ord=norm_ord)**norm_ord / N)

xx = range(1, max_d+1)
plt.plot(xx, err)
plt.plot(xx, test_err)
plt.plot([xx[0], xx[-1]], [lserr, lserr])
plt.plot([xx[0], xx[-1]], [test_lserr, test_lserr])
plt.xlabel('d')
plt.legend(('Training error', 'Test error', 'LS', 'LS test'))
np.argmax(test_err)
print(lsrank, np.max(ranks))
plt.plot(xx, ranks)
plt.plot((xx[0], xx[-1]),(lsrank, lsrank))
plt.legend(('rank after rf', 'original rank'))
plt.xlabel('d')<jupyter_output>99 100<jupyter_text>Playing with perturbed matrices<jupyter_code>n = 2000
m = 1000
r = 1
A = gen_rank(n,m,r)
std = 1

B = A + np.random.normal(size=(n, m), scale=std)
print(np.linalg.matrix_rank(B))

C = A + np.random.normal(size=m, scale=std)
print(np.linalg.matrix_rank(C))


print('RFN:')
D = A
rank = np.linalg.matrix_rank(D)
ranks = [rank]
print('%d/%d' % (rank, min(m, n)))
max_layers = 100
activation = np.cos
#activation = lambda X: np.maximum(X, 0)
#activation = lambda X: X**2 / np.linalg.norm(X**2, 2)
#activation = lambda X: 0.01 * X**2
while(np.linalg.matrix_rank(D) < min(n, m) and len(ranks) < max_layers):
    W = np.random.normal(size=(n,n), scale=std)
    b = np.random.normal(size=m, scale=std)
    D = np.matmul(W, D) + b
    D = activation(D)
    rank = np.linalg.matrix_rank(D)
    ranks.append(rank)
    print('%d/%d' % (rank, min(m, n)))
plt.plot(range(len(ranks)), ranks)<jupyter_output>1000
2
RFN:
1/1000
36/1000
1000/1000<jupyter_text>Fourier random features  Gaussian kernel<jupyter_code>#This represents functions in the RKHS of the Gaussian kernel k(x,x') = exp(-gamma||x-x'||^2) 
def make_fourier_rfm(d, m):
    W = np.random.normal(size=(m, d))
    b = np.random.uniform(low=0, high=2*math.pi, size=d)
    return lambda X: np.cos(np.matmul(X, W) + b)
def ls(A, b, reg=0):
    return np.linalg.solve(np.matmul(A.T, A) + reg * np.eye(A.shape[-1]), np.matmul(A.T, b))

def rf_fit(X, y, rfm, reg=0):
    Phi = rfm(X) #Nxd
    theta = ls(Phi, y, reg)
    yhat = np.matmul(Phi, theta)
    return yhat, theta, lambda Z: np.matmul(rfm(np.expand_dims(Z, 1)), theta)<jupyter_output><empty_output><jupyter_text>The Gaussian kernel is universal, which basically means we can represent all continuous functions (on compact sets)Let's try with polynomials (of scalar input)<jupyter_code>from numpy.polynomial.polynomial import Polynomial as Poly

domain = [-1, 1]

def make_poly(order=5):
    coeff = 2 * np.random.uniform(size=order+1) - 1
    poly = Poly(coeff, domain)
    return lambda x: poly(x).squeeze()
f = make_poly()
xx = np.linspace(domain[0], domain[1], num=1000)
plt.plot(xx, f(xx))
N = 100
noise = 0.1

X = 2 * np.random.uniform(size=(N, 1)) - 1
y = f(X) + noise * np.random.normal(size=N)
dmax = 20
reg = 0.1
results = []
ranks = []
for d in range(1, dmax+1):
    rfm = make_fourier_rfm(d, 1)
    results.append(rf_fit(X, y, rfm, reg))
    ranks.append(np.linalg.matrix_rank(rfm(X)))
xx = np.linspace(domain[0], domain[1], num=1000)

def update(d):
    plt.scatter(X.squeeze(), y)
    plt.plot(xx, results[d-1][2](xx), 'orange')
    axes = plt.gca()
    axes.set_xlim(domain)
    axes.set_ylim(min(y) - 0.1, max(y) + 0.1)
    
slider = widgets.interact(update, d=widgets.IntSlider(min=1, max=dmax, step=1, value=1))<jupyter_output><empty_output><jupyter_text>Rank of design matrix as a function of dimension<jupyter_code>plt.plot(range(1, dmax+1), ranks)<jupyter_output><empty_output><jupyter_text>Extreme machines<jupyter_code>def make_fourier_network(widths, input_dim):
    layers = []
    d_prev = input_dim
    for d in widths:
        W = np.random.normal(size=(d_prev, d))
        b = np.random.uniform(low=0, high=2*math.pi, size=d)
        layers.append({'W': W, 'b': b})
        d_prev = d
    def fourier_network(X):
        P = X
        for lay in layers:
            P = np.cos(np.matmul(P, lay['W']) + lay['b'])
        return P
    return fourier_network
layers = 2
dmax = 20
reg = 0.1
results = []
ranks = []
for d in range(1, dmax+1):
    rfm = make_fourier_network([d]*layers, 1)
    results.append(rf_fit(X, y, rfm, reg))
    ranks.append(np.linalg.matrix_rank(rfm(X)))
xx = np.linspace(domain[0], domain[1], num=1000)

def update(d):
    plt.scatter(X.squeeze(), y)
    plt.plot(xx, results[d-1][2](xx), 'orange')
    axes = plt.gca()
    axes.set_xlim(domain)
    axes.set_ylim(min(y) - 0.1, max(y) + 0.1)
    
slider = widgets.interact(update, d=widgets.IntSlider(min=1, max=dmax, step=1, value=1))<jupyter_output><empty_output><jupyter_text>Rank of design matrix as a function of width<jupyter_code>plt.plot(range(1, dmax+1), ranks)<jupyter_output><empty_output><jupyter_text>Gaussian RF for the Linear Rank-Deficient Scenario<jupyter_code>n = 2000
d = 100
r = 1
noise = 0.1
reg = 0.1
def gen_data(n, d, r):
    assert r <= d
    feats = np.random.normal(size=(n, d))
    proj = np.matmul(np.random.normal(size=(d,r)), np.random.normal(size=(r,d)))
    feats = np.matmul(feats, proj)
    feats = feats / np.max(feats)
    return feats
param = 2 * np.random.uniform(size=d) - 1
X = gen_data(n, d, r)
rank = np.linalg.matrix_rank(X)
y = np.matmul(X, param) + noise * np.random.normal(size=n)
param_ls = ls(X, y, reg)
paramerr = np.linalg.norm(param_ls - param)
yhat = np.matmul(X, param_ls)
mse_train = np.linalg.norm(yhat - y, ord=2)**2 / n
X_test = gen_data(n, d, r)
y_test = np.matmul(X_test, param) + noise * np.random.normal(size=n)
yhat_test = np.matmul(X_test, param_ls)
mse_test = np.linalg.norm(yhat_test - y_test, ord=2)**2 / n
rank, mse_train, mse_test, paramerr
rfn = make_fourier_network([100], 100)

yhat, param_rf, _ = rf_fit(X, y, rfn, reg)
mse_train = np.linalg.norm(yhat - y, ord=2)**2 / n
X_test = gen_data(n, d, r)
y_test = np.matmul(X_test, param) + noise * np.random.normal(size=n)
yhat_test = np.matmul(rfn(X_test), param_rf)
mse_test = np.linalg.norm(yhat_test - y_test, ord=2)**2 / n
rank = np.linalg.matrix_rank(rfn(X))
paramerr = None
rank, mse_train, mse_test, paramerr<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1778280190223747, 'fraction_numerical': 0.08762766040383566, 'mean_word_length': 4.2061688311688314, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '489218', 'license': 'MIT', 'n_tokens_mistral': 6578, 'n_tokens_neox': 5775, 'n_words': 1140, 'path': 'notebooks/Random Features.ipynb', 'repo_name': 'T3p/hidden-features', 'size': '12827'}