test.py

import numpy as np
import math
from scipy.spatial import distance

n = 2000
p = 800
r = 1000
print("r/n:", r/n)
print("p/n:", p/n)
print("theory result:", math.log(1+(n-p)/(r-p)))

TRUE_W = np.random.uniform(size=p)
TRUE_b = np.full(n, np.random.uniform(size=1))

X = np.random.uniform(size=(n, p))
noise = np.random.uniform(size=(n))
y = np.dot(X, TRUE_W) + TRUE_b + noise

r1 = np.linalg.inv(np.dot(X.T, X)).dot(X.T).dot(y)
#print(distance.cosine(r1, TRUE_W))

total_ve = 0.0
log_ve = 0.0
for i in range(10):
    S = np.random.normal(n, r)
    S = S / math.sqrt(r)
    X_S = np.dot(X.T, S).T
    y_s = np.dot(y.T, S).T
    r2 = np.linalg.inv(np.dot(X_S.T, X_S)).dot(X_S.T).dot(y_s)
#    print(distance.cosine(r2, TRUE_W))
    VE = np.linalg.norm(r2 - TRUE_W, ord=2)/np.linalg.norm(r1 - TRUE_W, ord=2)
    total_ve += VE
    log_ve += math.log(VE)
print(math.log(total_ve/10))
print(log_ve/10)