yaroslaw2.py

import numpy as np
import matplotlib.pyplot as plt

# Parameters
d = 5
n = 5
maxP = 10
n_samples = 10000

def formula_p1(n, d):
    return d*(1-1/d)**n
def formula_p2(n, d):
    return 1/3*d*(-((-1+d)*((-2+d)/d)**n)+(2+d)*((-1+d**2)/(d*(2+d)))**n)

def formula_p2_v2(n, d):
    mat = np.array([
        [ 1-2/d+2/(d*(d+2)), 1/(d*(d+2)) ],
        [ 2/(d*(d+2)), 1-2/d+1/(d*(d+2))]
    ])
    val = np.array([d, d**2])
    for _ in range(n):
        val = mat @ val
    return val[0]

def formula_p2_v3(n, d):
    mat = np.array([
        [ 1-2/d+1/(d*(d+2)), 1/(d*(d+2)), 1/(d*(d+2)) ],
        [ 1/(d*(d+2)), 1-2/d+1/(d*(d+2)), 1/(d*(d+2))],
        [0, 0, 1-1/d]
    ])
    val = np.array([d, d**2, d])
    for _ in range(n):
        val = mat @ val
    return val[0]


# Identity matrix
ii = np.eye(d)

# Initialize accumulators for averages
traces1_acc = np.zeros(maxP)
traces2_acc = np.zeros(maxP)
singular_sum_acc = np.zeros(maxP)

# Perform averaging over multiple random samples
for _ in range(n_samples):
    # Generate isotropic normalized random vectors
    isotropic = np.array([v / np.linalg.norm(v) for v in np.random.normal(size=(n, d))])

    # Compute A
    A_matrices = [ii - np.outer(v, v) for v in isotropic]
    A = np.linalg.multi_dot(A_matrices)
    B = A @ A.T

    # Compute powers of A
    A_powers = [A]
    B_powers = [B]
    for _ in range(1, maxP):
        A_powers.append(A_powers[-1] @ A)
        B_powers.append(B_powers[-1] @ B)

    # Accumulate trace values
    traces1_acc += np.array([np.trace(p) for p in A_powers])
    traces2_acc += np.array([np.trace(p) for p in B_powers])

    # Accumulate singular value calculations
    singular_values = np.linalg.svd(A, compute_uv=False)
    singular_sum_acc += np.array([np.sum(singular_values**(2 * p)) for p in range(1, maxP + 1)])

# Compute averages
traces1_avg = traces1_acc / n_samples
traces2_avg = traces2_acc / n_samples
singular_sum_avg = singular_sum_acc / n_samples

print(traces1_avg[:2])
print(traces2_avg[:2])
print(singular_sum_avg[:2])
print(formula_p1(n, d))
print(formula_p2(n, d))
print(formula_p2_v2(n, d))
print(formula_p2_v3(n, d))

# Trace plot
plt.plot(range(1, maxP + 1), traces1_avg, label="Avg Tr(A^p)", marker="o")
plt.plot(range(1, maxP + 1), traces2_avg, label="Avg Tr((AA^T)^p)", color="r", marker="+")
plt.plot(range(1, maxP + 1), singular_sum_avg, label="Avg Sum of Singular Values^(2p)", color="b")
plt.fill_between(range(1, maxP + 1), 0, traces1_avg, alpha=0.3)
plt.xlabel("p")
plt.ylabel("Average Trace")
plt.title(f"Average over {n_samples} samples, {d=}")
plt.legend()
plt.show()