import numpy as np

def freivalds_check(A, B, C, k=2):
    """
    Проверяет равенство A @ B = C
    с помощью алгоритма Фрейвалдса k раз.
    Возвращает True (скорее всего) или False (точно).
    """
    n = A.shape[0]
    for _ in range(k):
        # Случайный вектор (элементы -1 или 1)
        r = np.random.randint(0, 2, size=(n,)) * 2 - 1
        
        Br = B @ r      # O(n^2)
        C_r = C @ r     # O(n^2)
        ABr = A @ Br    # O(n^2)
        
        if not np.allclose(ABr, C_r, atol=1e-8):
            return False  # Точно не равны
    return True           # Скорее всего равны

# Демонстрация
n = 500
A = np.random.randn(n, n)
B = np.random.randn(n, n)
C = A @ B  # Истинное произведение

print("Проверяем истинное равенство A*B и C:")
print(freivalds_check(A, B, C, k=3))

# Немного подпортим C
C_perturbed = C.copy()
C_perturbed[0,0] += 1e-1

print("Проверяем заведомо неверное равенство:")
print(freivalds_check(A, B, C_perturbed, k=3))

Проверяем истинное равенство A*B и C:
True
Проверяем заведомо неверное равенство:
False

Рандомизированное матричное умножение

import numpy as np

def randomized_matrix_multiplication(A, B, k):
    """Приближенное матричное умножение AB с сэмплированием."""
    m, n = A.shape
    n_b, p = B.shape
    if n != n_b:
        raise ValueError("Несовместимые размеры матриц")

    # Вычисляем нормы для вероятностей (можно упростить до p_i = 1/n)
    col_norms_A = np.linalg.norm(A, axis=0)
    row_norms_B = np.linalg.norm(B, axis=1)
    probs = col_norms_A * row_norms_B
    sum_probs = np.sum(probs)
    if sum_probs == 0: # Если одна из матриц нулевая
         return np.zeros((m, p))
    probs /= sum_probs

    # Можно использовать равномерное распределение p_i = 1/n
    # probs = np.ones(n) / n

    C_approx = np.zeros((m, p))
    indices = np.random.choice(n, size=k, p=probs)
    # Если используем равномерное распределение:
    # indices = np.random.choice(n, size=k)
    # probs_selected = np.ones(k) / n

    probs_selected = probs[indices]

    for t in range(k):
        idx = indices[t]
        p_i = probs_selected[t]
        if p_i > 1e-9: # Избегаем деления на ноль
            A_col = A[:, idx:idx+1] # A^{(i_t)} как столбец (m x 1)
            B_row = B[idx:idx+1, :] # B_{(i_t)} как строка (1 x p)
            C_approx += (1.0 / (k * p_i)) * (A_col @ B_row)

    return C_approx

# Пример использования
m, n, p = 50, 100, 60
A = np.random.rand(m, n)
B = np.random.rand(n, p)
k = 500 # Количество сэмплов

C_exact = A @ B
C_approx = randomized_matrix_multiplication(A, B, k)

error = np.linalg.norm(C_exact - C_approx, 'fro') / np.linalg.norm(C_exact, 'fro')
print(f"Относительная ошибка Фробениуса: {error:.4f}")

Относительная ошибка Фробениуса: 0.0377

import numpy as np

def randomized_matmul(A, B, k):
    """
    Приближённое умножение матриц A (m x p) и B (p x n).
    Возвращает матрицу C ~ A @ B.
    """
    m, p = A.shape
    p2, n = B.shape
    assert p == p2
    
    # Вероятности берём пропорционально нормам столбцов A * строк B
    col_norms = np.linalg.norm(A, axis=0)  # размер p
    row_norms = np.linalg.norm(B, axis=1)  # размер p
    weights = col_norms * row_norms
    weights = weights / np.sum(weights)
    
    # Сэмплим k индексов без возвращения
    idx = np.random.choice(p, size=k, replace=False, p=weights)
    
    # Формируем сумму
    C_approx = np.zeros((m, n))
    for i in idx:
        # A^{(i)} -- i-ый столбец A, B_{(i)} -- i-ая строка B
        # Учитываем масштабирование 1/(k * p_i)
        C_approx += (1.0 / (k * weights[i])) * np.outer(A[:, i], B[i, :])
    
    return C_approx

# Тест
m, p, n = 100, 500, 100
A = np.random.randn(m, p)
B = np.random.randn(p, n)
C_true = A @ B

# k намного меньше p
k = 500
C_approx = randomized_matmul(A, B, k)

error_rel = np.linalg.norm(C_true - C_approx, 'fro') / np.linalg.norm(C_true, 'fro')
print(f"Относительная ошибка (F-норма): {error_rel:.4f}")

Относительная ошибка (F-норма): 0.0902

Оценка следа матрицы

import numpy as np

def hutchinson_trace_estimator(A, k):
    """Оценка следа Tr(A) методом Хатчинсона."""
    n = A.shape[0]
    trace_estimates = []
    for _ in range(k):
        w = np.random.choice([-1.0, 1.0], size=(n, 1))
        trace_estimate = (w.T @ A @ w).item() # w^T A w
        trace_estimates.append(trace_estimate)
    return np.mean(trace_estimates)

def girard_trace_estimator(A, k):
    """Оценка следа Tr(A) методом Жирарда."""
    n = A.shape[0]
    trace_estimates = []
    for _ in range(k):
        w = np.random.randn(n, 1) # w ~ N(0, I)
        trace_estimate = (w.T @ A @ w).item() # w^T A w
        trace_estimates.append(trace_estimate)
    return np.mean(trace_estimates)

# Пример использования
n = 100
# Создаем симметричную положительно определенную матрицу
A_sym = np.random.rand(n, n)
A_sym = 0.5 * (A_sym + A_sym.T)
A_sym = A_sym + n * np.eye(n) # Гарантируем положительную определенность

exact_trace = np.trace(A_sym)
k = 100 # Количество векторов для сэмплирования

hutchinson_est = hutchinson_trace_estimator(A_sym, k)
girard_est = girard_trace_estimator(A_sym, k)

print(f"Точный след: {exact_trace:.4f}")
print(f"Оценка Хатчинсона (k={k}): {hutchinson_est:.4f}, Ошибка: {abs(hutchinson_est - exact_trace):.4f}")
print(f"Оценка Жирарда (k={k}): {girard_est:.4f}, Ошибка: {abs(girard_est - exact_trace):.4f}")

# Оценка intdim
norm_f_sq = np.linalg.norm(A_sym, 'fro')**2
intdim = (exact_trace**2) / norm_f_sq
print(f"Внутренняя размерность (intdim): {intdim:.4f} (макс. {n})")

Точный след: 10055.8264
Оценка Хатчинсона (k=100): 10048.4577, Ошибка: 7.3687
Оценка Жирарда (k=100): 10072.2824, Ошибка: 16.4561
Внутренняя размерность (intdim): 99.7131 (макс. 100)

import numpy as np

def hutchinson_trace_estimator(A_mv, n, k=10):
    """
    Оценивает Tr(A), не зная A явно, 
    но имея функцию A_mv(x), вычисляющую A@x.
    
    Параметры:
    -----------
    A_mv : function
        Функция, принимающая на вход x и возвращающая A@x.
    n : int
        Размерность матрицы A (n x n).
    k : int
        Число случайных векторов для усреднения.
    """
    estimates = []
    for _ in range(k):
        # Случайный радемахеровский вектор
        r = np.random.randint(0, 2, size=n)*2 - 1
        Ar = A_mv(r)
        estimates.append(r @ Ar)  # скаляр w^T (A w)
    return np.mean(estimates)

# Пример: матрица A и явный A_mv
n = 30
A = np.random.randn(n, n)
A = 0.5*(A + A.T)  # Сделаем её симметричной для наглядности
true_trace = np.trace(A)

# Создаём функцию A_mv
def A_mv(x):
    return A @ x

k = 8000
trace_hat = hutchinson_trace_estimator(A_mv, n, k)
print(f"Истинный след:      {true_trace:.2f}")
print(f"Оценка Хатчинсона: {trace_hat:.2f}")

Истинный след:      -9.50
Оценка Хатчинсона: -9.18

import numpy as np

def girard_trace_estimator(A_mv, n, k=10):
    """
    Оценивает Tr(A), используя гауссовские случайные векторы.
    """
    estimates = []
    for _ in range(k):
        w = np.random.randn(n)  # гауссовский вектор
        Aw = A_mv(w)
        estimates.append(w @ Aw)
    return np.mean(estimates)

def intrinsic_dimension(A):
    # Предполагаем, что A задана явно, чтобы просто продемонстрировать
    trace_val = np.trace(A)
    fro_norm = np.linalg.norm(A, 'fro')
    return trace_val / fro_norm

# Демонстрация
n = 200
A = np.random.randn(n, n)
A = A @ A.T  # делаем A = A^T A, чтобы получилась ПД матрица
true_trace = np.trace(A)

def A_mv(x):
    return A @ x

k = 30
trace_est_girard = girard_trace_estimator(A_mv, n, k)
idim = intrinsic_dimension(A)

print(f"Истинный след:          {true_trace:.2f}")
print(f"Оценка Жирара (k={k}):  {trace_est_girard:.2f}")
print(f"Intrinsic dimension:    {idim:.2f}")

Истинный след:          39814.18
Оценка Жирара (k=30):  40171.33
Intrinsic dimension:    9.99

Рандомизированный SVD

import numpy as np

def randomized_svd(A, rank, oversampling=10):
    """Вычисляет приближенное SVD ранга 'rank'."""
    m, n = A.shape
    ell = rank + oversampling # Целевой ранг + доп. сэмплирование

    # Шаг 1: Генерируем случайную матрицу
    Omega = np.random.randn(n, ell)

    # Шаг 2: Вычисляем Y = A * Omega
    Y = A @ Omega

    # Шаг 3: Находим ортонормированный базис Q для Y (QR разложение)
    Q, _ = np.linalg.qr(Y) # Q имеет размер (m x ell)

    # Шаг 4: Проецируем A на базис Q: B = Q^T * A
    B = Q.T @ A # B имеет размер (ell x n)

    # Шаг 5: Вычисляем SVD для маленькой матрицы B
    U_tilde, s, Vh = np.linalg.svd(B, full_matrices=False) # U_tilde (ell x ell), s (ell,), Vh (ell x n)

    # Шаг 6: Формируем U для исходной матрицы A
    U = Q @ U_tilde # U имеет размер (m x ell)

    # Шаг 7: Обрезаем до нужного ранга k
    U_k = U[:, :rank]
    s_k = s[:rank]
    Vh_k = Vh[:rank, :]

    return U_k, s_k, Vh_k

# Пример использования
m, n = 1000, 500
A = np.random.rand(m, n) @ np.random.rand(n, n) # Создаем матрицу с падающими сингулярными числами

target_rank = 10

U_r, s_r, Vh_r = randomized_svd(A, target_rank)
A_approx_r = U_r @ np.diag(s_r) @ Vh_r

# Для сравнения - точное SVD
U_exact, s_exact, Vh_exact = np.linalg.svd(A, full_matrices=False)
A_approx_exact_k = U_exact[:, :target_rank] @ np.diag(s_exact[:target_rank]) @ Vh_exact[:target_rank, :]

error_rand = np.linalg.norm(A - A_approx_r, 'fro') / np.linalg.norm(A, 'fro')
error_exact_k = np.linalg.norm(A - A_approx_exact_k, 'fro') / np.linalg.norm(A, 'fro')

print(f"Ранг приближения: {target_rank}")
print(f"Ошибка рандомизированного SVD: {error_rand:.4e}")
print(f"Ошибка точного SVD ранга k:     {error_exact_k:.4e}")
print(f"Первые 5 синг. чисел (рандом): {s_r[:5]}")
print(f"Первые 5 синг. чисел (точно):  {s_exact[:5]}")

Ранг приближения: 10
Ошибка рандомизированного SVD: 1.4632e-02
Ошибка точного SVD ранга k:     1.4126e-02
Первые 5 синг. чисел (рандом): [88663.64610265   107.31285841   104.0445584    102.02981015
   101.02648135]
Первые 5 синг. чисел (точно):  [88663.91400357   138.4687148    135.01731642   134.03101164
   132.74498632]

import numpy as np

def randomized_svd(A, k, p=5):
    """
    Возвращает приближение SVD матрицы A: A ~ U @ np.diag(S) @ V.
    k - желаемый ранг, p - оверсэмплинг.
    """
    m, n = A.shape
    # 1. Случайная матрица G
    G = np.random.randn(n, k + p)
    # 2. Y = A G
    Y = A @ G  # размер (m x (k+p))
    # 3. QR-разложение
    Q, _ = np.linalg.qr(Y, mode='reduced')  # Q: (m x (k+p))
    # 4. B = Q^T A
    B = Q.T @ A  # размер ((k+p) x n)
    # 5. Точное SVD матрицы B
    u_hat, s, v = np.linalg.svd(B, full_matrices=False)
    # 6. U = Q u_hat
    U = Q @ u_hat  # размер (m x (k+p))
    # Обрежем до ранга k
    return U[:, :k], s[:k], v[:k, :]

# Демонстрация
m, n, k = 1000, 300, 50
A = np.random.randn(m, k) @ np.random.randn(k, n)  # матрица ранга k (точно)
A += 0.01*np.random.randn(m, n)  # чуть портим

U_approx, S_approx, V_approx = randomized_svd(A, k=50, p=30)
A_approx = U_approx @ np.diag(S_approx) @ V_approx

error_fro = np.linalg.norm(A - A_approx, 'fro') / np.linalg.norm(A, 'fro')
print(f"Относительная ошибка (F-регион): {error_fro:.4e}")

Относительная ошибка (F-регион): 2.0224e-03

Kaczmarz algorithm

import numpy as np

def randomized_kaczmarz(A, b, iters=1000):
    """
    Решает систему A x = b методом Kaczmarz.
    Возвращает приблизительное решение x.
    """
    m, n = A.shape
    x = np.zeros(n)
    row_norms = np.sum(A*A, axis=1)
    row_norms_sum = np.sum(row_norms)
    
    for _ in range(iters):
        # Выбираем строку i случайно, пропорц. ее норме
        i = np.random.choice(m, p=row_norms/row_norms_sum)
        a_i = A[i, :]
        residual = a_i.dot(x) - b[i]
        x = x - (residual / (a_i.dot(a_i))) * a_i
    return x

# Демонстрация
m, n = 500, 300
A = np.random.randn(m, n)
x_true = np.random.randn(n)
b = A @ x_true

x_est = randomized_kaczmarz(A, b, iters=5000)
error = np.linalg.norm(x_est - x_true)/np.linalg.norm(x_true)
print(f"Относительная ошибка в решении: {error:.4e}")

Относительная ошибка в решении: 7.3028e-02