singulars.py

import numpy as np
import matplotlib.pyplot as plt

# Parameters
d = 10
maxP = 10
n_samples = 100

# 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=(d, 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

# 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()