<jupyter_start><jupyter_text>Generate ill-conditioned matrix, ```eps_cond``` controls condition number<jupyter_code>m = 10000
n = 100
eps_cond = 1e-6
seed = 42
np.random.seed(seed)

A = np.random.randn(m, n)
U, s, V = np.linalg.svd(A, full_matrices=False)
A = U.dot(np.diag(np.maximum(s - np.min(s) + eps_cond, 0))).dot(V)
s = np.linalg.svd(A, compute_uv=False, full_matrices=False)
print("Cond number = {}".format(max(s) / min(s)))
kappa = np.linalg.norm(A, "fro") / min(s)
print("Convergence factor for a consistent linear system = {}".format(kappa))
x_true = np.random.randn(n)<jupyter_output>Cond number = 19605065.815627676
Convergence factor for a consistent linear system = 112609049.85223769<jupyter_text>Construct consistent or noisy linear system, control it with a noise level ```sigma```<jupyter_code>sigma = 0.
b = A.dot(x_true) + sigma * np.random.randn(m)
b = b[:, np.newaxis]
print(A.shape, b.shape)<jupyter_output>(10000, 100) (10000, 1)<jupyter_text>Construct preconditioners from sketched matrix<jupyter_code>gamma_range = list(range(1, 5))

preprocess_time = { "PrecRKM, $\gamma = {}$".format(gamma): 0 for gamma in gamma_range}

sketch_sizes = [gamma * n for gamma in gamma_range]
P = {name: [] for name in preprocess_time}
# From some random submatrix
P = {}
for k, gamma in enumerate(gamma_range):
    print("Sketch size = {}".format(n * gamma))
    st = time.time()
    sketch_idx = np.sort(np.random.choice(np.arange(m), sketch_sizes[k], replace=False))
    Q, R_qr = np.linalg.qr(A[sketch_idx, :])
    P["PrecRKM, $\gamma = {}$".format(gamma)] = np.linalg.inv(R_qr)
    preprocess_time["PrecRKM, $\gamma = {}$".format(gamma)] = time.time() - st<jupyter_output>Sketch size = 100
Sketch size = 200
Sketch size = 300
Sketch size = 400<jupyter_text>Defines parameters of the compared methods<jupyter_code>max_iter = 300000
tol = 1e-6
x0 = np.zeros(n)
# Maximum number of seconds to run the method
max_time = 1.5
# How frequent the current iterand is stored
log_interval = 10<jupyter_output><empty_output><jupyter_text>Compare preconditioned Kaczmarz method with no preconditioning<jupyter_code>conv_x_res = {"RKM": []}
conv_time_res = {"RKM": []}
res = methods.classical_kaczmarz(A, b, x0, max_iter, max_time, log_interval)
conv_x_res["RKM"] = res["conv_x"]
conv_time_res["RKM"] = res["conv_time"]
for name in preprocess_time:
    print(name)
    res = methods.precond_kaczmarz(A, b, P[name], 
                             x0, max_iter, max_time, 
                             log_interval, preprocess_time[name])
    conv_x_res[name] = res["conv_x"]
    conv_time_res[name] = res["conv_time"]
res = methods.extended_kaczmarz(A, b, x0, max_iter, max_time, log_interval)
conv_x_res["REK"] = res["conv_x"]
conv_time_res["REK"] = res["conv_time"]
res = methods.rdbk(A, b, x0, 10, 10, max_iter, max_time, log_interval)
conv_x_res["Block REK"] = res["conv_x"]
conv_time_res["Block REK"] = res["conv_time"]
relative_residual_norm = {}
for name in conv_x_res:
    relative_residual_norm[name] = []
    if P.get(name) is None:
        for x in conv_x_res[name]:
            relative_residual_norm[name].append(np.linalg.norm(A @ x - b) / np.linalg.norm(b))
    else:
        for x in conv_x_res[name]:
            relative_residual_norm[name].append(np.linalg.norm(A @ (P[name] @ x) - b) / np.linalg.norm(b))<jupyter_output><empty_output><jupyter_text>Plots Relative error vs. number of iterations<jupyter_code>plt.figure(figsize=(10, 8))

for name in relative_residual_norm:
    iter_range = [j * log_interval for j in range(len(conv_x_res[name]))]
    plt.semilogy(iter_range, relative_residual_norm[name], label=name, linewidth=3)
plt.ylabel(r"$\frac{\|A\hat{P}_Rx_k - b\|_2}{\|b\|_2}$", fontsize=36)
plt.xlabel(r"Number of iterations, $k$", fontsize=24)
plt.legend(fontsize=18)
plt.grid(True)
_ = plt.yticks(fontsize=18)
_ = plt.xticks(fontsize=18)<jupyter_output><empty_output><jupyter_text>Relative error vs. running time + time to construct preconditioners<jupyter_code>plt.figure(figsize=(10, 8))

for name in relative_residual_norm:
    print(name)
    if conv_time_res[name][-1] == conv_time_res[name][::log_interval][-1]:
        time_range = conv_time_res[name][::log_interval]
    else:
        time_range = conv_time_res[name][::log_interval] + [conv_time_res[name][-1]]
    plt.semilogy(time_range, relative_residual_norm[name], label=name, linewidth=3)

plt.ylabel(r"$\frac{\|A\hat{P}_Rx_k - b\|_2}{\|b\|_2}$", fontsize=34)
plt.xlabel(r"CPU time, sec", fontsize=34)
plt.legend(fontsize=20)
_ = plt.yticks(fontsize=28)
_ = plt.xticks(fontsize=28)
plt.grid(True)<jupyter_output>RKM
PrecRKM, $\gamma = 1$
PrecRKM, $\gamma = 2$
PrecRKM, $\gamma = 3$
PrecRKM, $\gamma = 4$
REK
Block REK
{'fraction_non_alphanumeric': 0.1547442156654638, 'fraction_numerical': 0.02908087454892804, 'mean_word_length': 5.056555269922879, 'pattern_counts': {'":': 2, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '477143', 'license': 'MIT', 'n_tokens_mistral': 1996, 'n_tokens_neox': 1834, 'n_words': 364, 'path': 'examples/random_data.ipynb', 'repo_name': 'amkatrutsa/preckacz', 'size': '4711'}