Я реализовал так называемый countSketch в Python (стр. 17: https://arxiv.org/pdf/1411.4357.pdf), но моей реализации в настоящее время не хватает производительности. Алгоритм состоит в том, чтобы вычислить продукт SA, где A является n x d матрица, S - это m x n матрица, определенная следующим образом: для каждого столбца S равномерно случайным образом выберите строку (хэш-корзину) из строк m и для этой данной строки, равномерно при случайном выборе +1 или -1, поэтому S - это матрица с ровно одним ненулевым значением в каждом столбце, а в остальном все нулями.
Мое намерение состоит в том, чтобы вычислить SA в потоковом режиме , читая записииз A. Идея для моей реализации заключается в следующем: наблюдаем последовательность троек (i,j,A_ij) и возвращаем последовательность (h(i), j, s(i)A_ij), где: - h(i) - это хэш-контейнер (строка матрицы, выбранная случайным образом из m возможных строк S - s(i) - это функция случайного знака, как описано выше. Я предположил, что матрица находится в порядке строк , так что первая строка AПолностью прибывает до того, как прибудет следующая строка A, потому что это ограничивает количество вызовов, которые мне нужны для выбора случайного сегмента, или необходимость использовать хеш-библиотеку.Я также предположил, что число ненулевых записей (или длина входного потока) известно, так что я могу эффективно кодировать итерацию.
Моя проблема в том, что матрица должна вычислять (1+error)*||Ax||^2 <= ||SAx||^2 <= (1+error)*||Ax||^2, а также иметьразница в нормах Фробениуса между A^T S^T S A и A^T A невелика.Однако, хотя моя реализация первого условия кажется верной, последнее всегда слишком мало.Мне было интересно, есть ли очевидная причина для этого, которую я пропускаю, потому что кажется, что она недооценивает последнее количество.
Я открыт для обратной связи по изменению кода, если есть очевидные улучшения, которые необходимо сделать.Один вызов np.choice сделан для того, чтобы избавить от необходимости просматривать (потенциально большой) массив или хеш-таблицу, содержащую хеш-строки для каждой строки, и поскольку матрица находится в порядке строк, мы можем просто сохранить хеш-код для этой строки доновая строка видна.
nb.Если вы не хотите запускать с помощью numba, просто закомментируйте импорт и декоратор функций, и он будет работать в стандартном режиме numpy / scipy.
import numpy as np
import numpy.random as npr
import scipy.sparse as sparse
from scipy.sparse import coo_matrix
import numba
from numba import jit
@jit(nopython=True) # comment this if want just numpy
def countSketch(input_rows, input_data,
input_nnz,
sketch_size, seed=None):
'''
input_rows: row indices for data (can be repeats)
input_data: values seen in row location,
input_nnz : number of nonzers in the data (can replace with
len(input_data) but avoided here for speed)
sketch_size: int
seed=None : random seed
'''
hashed_rows = np.empty(input_rows.shape,dtype=np.int32)
current_row = 0
hash_val = npr.choice(sketch_size)
sign_val = npr.choice(np.array([-1.0,1.0]))
#print(hash_val)
hashed_rows[0] = hash_val
#print(hash_val)
for idx in np.arange(input_nnz):
print(idx)
row_id = input_rows[idx]
data_val = input_data[idx]
if row_id == current_row:
hashed_rows[idx] = hash_val
input_data[idx] = sign_val*data_val
else:
# make new hashes
hash_val = npr.choice(sketch_size)
sign_val = npr.choice(np.array([-1.0,1.0]))
hashed_rows[idx] = hash_val
input_data[idx] = sign_val*data_val
return hashed_rows, input_data
def sort_row_order(input_data):
sorted_row_column = np.array((input_data.row,
input_data.col,
input_data.data))
idx = np.argsort(sorted_row_column[0])
sorted_rows = np.array(sorted_row_column[0,idx], dtype=np.int32)
sorted_cols = np.array(sorted_row_column[1,idx], dtype=np.int32)
sorted_data = np.array(sorted_row_column[2,idx], dtype=np.float64)
return sorted_rows, sorted_cols, sorted_data
if __name__=="__main__":
import time
from tabulate import tabulate
matrix = sparse.random(1000, 50, 0.1)
x = np.random.randn(matrix.shape[1])
true_norm = np.linalg.norm(matrix@x,ord=2)**2
tidy_data = sort_row_order(matrix)
sketch_size = 300
start = time.time()
hashed_rows, sketched_data = countSketch(tidy_data[0],\
tidy_data[2], matrix.nnz,sketch_size)
duration_slow = time.time() - start
S_A = sparse.coo_matrix((sketched_data, (hashed_rows,matrix.col)))
approx_norm_slow = np.linalg.norm(S_A@x,ord=2)**2
rel_error_slow = approx_norm_slow/true_norm
#print("Sketch time: {}".format(duration_slow))
start = time.time()
hashed_rows, sketched_data = countSketch(tidy_data[0],\
tidy_data[2], matrix.nnz,sketch_size)
duration = time.time() - start
#print("Sketch time: {}".format(duration))
S_A = sparse.coo_matrix((sketched_data, (hashed_rows,matrix.col)))
approx_norm = np.linalg.norm(S_A@x,ord=2)**2
rel_error = approx_norm/true_norm
#print("Relative norms: {}".format(approx_norm/true_norm))
print(tabulate([[duration_slow, rel_error_slow, 'Yes'],
[duration, rel_error, 'No']],
headers=['Sketch Time', 'Relative Error', 'Dry Run'],
tablefmt='orgtbl'))