#python #numpy #scipy
#python #numpy #scipy
Вопрос:
У меня есть два numpy-вектора a и b, и я хочу построить (идеально разреженную) матрицу C с помощью чего-то вроде:
C[i,j] = (a[j]>b[j])*(a[i] b[j])
Таким образом, только те элементы, в которых соблюдается определенное соотношение между j-индексами, будут вычисляться с использованием некоторых других элементов. Опять же, это легко сделать с помощью цикла i amp; j, но я хотел бы знать, есть ли более эффективный способ «numpy / scipy»? Проще говоря, я не мог понять, как справиться с различиями в используемых индексах.
Заранее спасибо!
Ответ №1:
Я должен был настаивать на [MCVE], но я сделаю один
Для массивов одинакового размера:
In [232]: a = np.random.randint(0,10,10)
In [233]: b = np.random.randint(0,10,10)
a[j]>b[j] часть вашего вычисления:
In [234]: a>b
Out[234]:
array([ True, False, False, False, False, True, False, False, False,
False])
часть суммирования, a[i] b[j] :
In [235]: a[:,None] b
Out[235]:
array([[ 9, 9, 10, 14, 13, 11, 11, 12, 9, 14],
[ 4, 4, 5, 9, 8, 6, 6, 7, 4, 9],
[ 7, 7, 8, 12, 11, 9, 9, 10, 7, 12],
[10, 10, 11, 15, 14, 12, 12, 13, 10, 15],
[12, 12, 13, 17, 16, 14, 14, 15, 12, 17],
[12, 12, 13, 17, 16, 14, 14, 15, 12, 17],
[ 6, 6, 7, 11, 10, 8, 8, 9, 6, 11],
[11, 11, 12, 16, 15, 13, 13, 14, 11, 16],
[ 7, 7, 8, 12, 11, 9, 9, 10, 7, 12],
[12, 12, 13, 17, 16, 14, 14, 15, 12, 17]])
Установка желаемых столбцов в 0:
In [236]: _[:,np.nonzero(a>b)] = 0
In [237]: _
Out[237]:
array([[ 0, 9, 10, 14, 13, 0, 11, 12, 9, 14],
[ 0, 4, 5, 9, 8, 0, 6, 7, 4, 9],
[ 0, 7, 8, 12, 11, 0, 9, 10, 7, 12],
[ 0, 10, 11, 15, 14, 0, 12, 13, 10, 15],
[ 0, 12, 13, 17, 16, 0, 14, 15, 12, 17],
[ 0, 12, 13, 17, 16, 0, 14, 15, 12, 17],
[ 0, 6, 7, 11, 10, 0, 8, 9, 6, 11],
[ 0, 11, 12, 16, 15, 0, 13, 14, 11, 16],
[ 0, 7, 8, 12, 11, 0, 9, 10, 7, 12],
[ 0, 12, 13, 17, 16, 0, 14, 15, 12, 17]])
Упс, я переключил это, я должен был установить для других столбцов значение 0.
Но нам не нужно делать это отдельными шагами:
In [238]: (a>b)*(a[:,None] b)
Out[238]:
array([[ 9, 0, 0, 0, 0, 11, 0, 0, 0, 0],
[ 4, 0, 0, 0, 0, 6, 0, 0, 0, 0],
[ 7, 0, 0, 0, 0, 9, 0, 0, 0, 0],
[10, 0, 0, 0, 0, 12, 0, 0, 0, 0],
[12, 0, 0, 0, 0, 14, 0, 0, 0, 0],
[12, 0, 0, 0, 0, 14, 0, 0, 0, 0],
[ 6, 0, 0, 0, 0, 8, 0, 0, 0, 0],
[11, 0, 0, 0, 0, 13, 0, 0, 0, 0],
[ 7, 0, 0, 0, 0, 9, 0, 0, 0, 0],
[12, 0, 0, 0, 0, 14, 0, 0, 0, 0]])
Это вычисление сильно зависит от broadcasting . Если вы не знакомы с этим, это не будет иметь особого смысла.
Массивы:
In [239]: a
Out[239]: array([5, 0, 3, 6, 8, 8, 2, 7, 3, 8])
In [240]: b
Out[240]: array([4, 4, 5, 9, 8, 6, 6, 7, 4, 9])
Тот факт, что многие (или только некоторые) столбцы будут равны 0, не имеет большого значения при выполнении numpy вычислений всего массива. Обычно проще и быстрее передавать целые массивы.
Вы можете выполнять итерации по столбцам и просто выполнять добавление для выбранных. Нет необходимости выполнять итерации по строкам. Но обычно это будет медленнее, чем [238]:
In [245]: c = np.zeros((10,10),int)
...: for j in range(10):
...: if a[j]>b[j]:
...: c[:,j] = a b[j]
или перемещение сравнения за пределы цикла:
In [249]: c = np.zeros((10,10),int)
...: m = a>b
...: for j,v in enumerate(m):
...: if v:
...: c[:,j] = a b[j]
Еще лучше:
In [253]: c = np.zeros((10,10),int)
...: idx = np.nonzero(a>b)[0]
...: print(idx)
...: c[:,idx] = a[:,None] b[idx]
Для некоторой степени разреженности это последнее может быть быстрее, но не здесь:
In [256]: %%timeit
...: c = np.zeros((10,10),int)
...: idx = np.nonzero(a>b)[0]
...: c[:,idx] = a[:,None] b[idx]
17.7 µs ± 772 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [258]: timeit (a>b)*(a[:,None] b)
10.7 µs ± 139 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Ответ №2:
Нет ничего плохого в использовании циклов. Есть одно свойство, которое вы можете использовать, подумав: Условие a[j]>b[j] сообщает нам, может ли весь столбец быть равен нулю, поэтому, если мы используем scipy.sparse.coo_matrix((data,(i,j))) , имеет смысл создать внешний цикл над столбцами:
data=[]
i=[]
j=[]
for jx = ...:
if a[j]>b[j]:
for ix = ...: # this whole loop can be skipped if the condition is not satisfied
data.append(a[i] a[j])
i.append(ix)
j.append(jx)
Комментарии:
1. Спасибо за предложение. Я должен был подчеркнуть, что этот процесс будет выполняться очень часто, и, насколько я думал, добавление является одним из довольно медленных процессов, не так ли? И это также будет выполняться однопоточно, верно?
2. В этом случае должно быть легко распараллелить цикл
jx. Добавление происходит не медленно, поскольку это операция O (1), но если это проблема, вы также можете рассмотреть возможность предварительного выделенияi,jкакnp.arrays и в зависимости от того, какaиbраспределены, вы могли бы ускорить его, векторизовав второй цикл, возможно, даже сделав уступку разреженности, которая для каждогоj, гдеa[j]>b[j]вы простовведите полный столбец. Но это то, что можете решить только вы, поскольку мы недостаточно знаем о вашей проблеме и ваших инструментах!3. В зависимости от вашего варианта использования может даже не потребоваться реализовать создание этих матриц, но, возможно, реализовать только метод, который вычисляет результат умножения матрицы. Но, как я уже сказал, это не то, о чем мы можем судить со стороны. В качестве альтернативы вы можете попытаться векторизовать все, что может быть даже быстрее, но это опять же зависит от того, насколько велики
aиb, какое у них распределение чисел и т. Д. То, что я предоставил здесь, было именно тем, что я бы счел наиболее «numpy / scipy»-способом.4. ОП сказал, что у него были
numpyмассивы. Обычно мы стараемся избегать циклов (python) при работе с массивами.5. @hpaulj В целом это верно, однако в данном случае это звучит так, как будто мы имеем дело с разреженными данными. Допустим, условие выполняется только в 0,01% случаев
a[j]>b[j], иaиbдостаточно велики. Тогда нам, возможно, действительно придется работать с циклами, поскольку вся матрица может даже не поместиться в память. И память, вероятно, является проблемой, если OP уже хочет работать с разреженными структурами данных. Но это всего лишь предположение, поскольку мы не знаем точной ситуации, в которой находится OP.