Numpy, эффективно вставить одну матрицу в другую?

Question

Я пытаюсь сделать внешнее произведение двух векторов более эффективным, удаляя нулевые элементы, выполняя внешнее произведение, а затем увеличивая полученную матрицу рядами нулей или вставляя в нулевую матрицу.(Разреживание матрицы с помощью scipy на самом деле не работает, так как стоимость конвертации высока, и я делаю это снова и снова.)

import numpy
dim = 100
vec = np.random.rand(1, dim)
mask = np.flatnonzero(vec > 0.8)
vec_sp = vec[:, mask]
mat_sp = vec_sp.T * vec_sp # This is faster than dot product
# Enlarge matrix or insert into zero matrix

Поскольку это внешнее произведение двух векторов, я знаю нольстроки и столбцы в исходной матрице, они являются индексами в переменной маски.Чтобы увидеть это,

a = np.array(((1,0,2,0))).reshape(1,-1)
a.T * a
>> array([[1, 0, 2, 0],
       [0, 0, 0, 0],
       [2, 0, 4, 0],
       [0, 0, 0, 0]])

Я попробовал два разных решения: одно с использованием метода insert numpy и добавление метода к переменной mat_sp.Все это становится циклом и действительно медленным.

for val in mask:
    if val < mat_sp.shape[0]:
        mat_sp = np.insert(mat_sp, val, values=0, axis=1)
        mat_sp = np.insert(mat_sp, val, values=0, axis=0)
    else:
        mat_sp = np.append(mat_sp, values=np.zeros((mat_sp.shape[0], 1)), axis=1)
        mat_sp = np.append(mat_sp, values=np.zeros((1, mat_sp.shape[1])), axis=0)

Другой подход заключается в создании нулевой матрицы размером dim x dim и последующем создании гигантского вектора индекса из маски через два цикла for.А затем с помощью вектора индекса вставить матричное умножение в нулевую матрицу.Тем не менее, это также очень медленно.

Любые идеи или идеи, которые могли бы эффективно решить проблему, были бы полезны, так как произведение разреженной матрицы занимает 2/3 времени не разреженного.

---

Используя пример@hjpaul мы получаем следующий код сравнения

import numpy as np
dims = 400

def test_non_sparse():
    vec = np.random.rand(1, dims)
    a = vec.T * vec

def test_sparse():  
    vec = np.random.rand(1, dims)
    idx = np.flatnonzero(vec>0.75)
    oprod = vec[:,idx].T * vec[:,idx]
    vec_oprod = np.zeros((dims, dims))
    vec_oprod[idx[:,None], idx] = oprod


if __name__ == '__main__':
    import timeit
    print('Non sparse:',timeit.timeit("test_non_sparse()", setup="from __main__ import test_non_sparse", number=10000))
    print('Sparse:',timeit.timeit("test_sparse()", setup="from __main__ import test_sparse", number=10000))

Код, конечно, дает улучшение в зависимости от размеров векторов и числа нулей.Свыше 300 измерений и около 70% нулей дают умеренное улучшение скорости, которое увеличивается с увеличением числа нулевых элементов и размеров.Если матрицы и маска одинаковы снова и снова, несомненно, можно добиться большего ускорения.

(Моя ошибка при выполнении логической индексации была idx вместо idx[:,None])

Answer 1

Самый быстрый способ вставить одну матрицу в другую с помощью индексов.

Для иллюстрации с вашим внешним продуктом:

In [94]: a = np.array(((1,0,2,0)))
In [95]: idx = np.where(a>0)[0]
In [96]: idx
Out[96]: array([0, 2])
In [97]: a1 = a[idx]

Внешний продукт массива конденсатов:

In [98]: a2 = a1[:,None]*a1
In [99]: a2
Out[99]: 
array([[1, 2],
       [2, 4]])

Настройте массив результатов и используйте блочную индексацию, чтобы определить, куда идут значения a2:

In [100]: res = np.zeros((4,4),int)
In [101]: res[idx[:,None], idx] = a2
In [102]: res
Out[102]: 
array([[1, 0, 2, 0],
       [0, 0, 0, 0],
       [2, 0, 4, 0],
       [0, 0, 0, 0]])

Прямой внешний вид неконденсированного массива:

In [103]: a[:,None]*a
Out[103]: 
array([[1, 0, 2, 0],
       [0, 0, 0, 0],
       [2, 0, 4, 0],
       [0, 0, 0, 0]])
In [104]: np.outer(a,a)
Out[104]: 
array([[1, 0, 2, 0],
       [0, 0, 0, 0],
       [2, 0, 4, 0],
       [0, 0, 0, 0]])

Если a было 2d, (n, 1), это внешнее можно записать как np.dot(a.T,a).dot включает в себя сумму, в данном случае измерения размера 1.

Я думаю, a пришлось бы довольно разреженным, чтобы извлечь выгоду из этой дополнительной работы по индексированию.Я обнаружил, что с разреженными разреженными матрицами разреженность порядка 1% дает преимущество в скорости, даже если матрицы предварительно изготовлены.

---

In [105]: from scipy import sparse
In [106]: A = sparse.csr_matrix(a)
In [107]: A
Out[107]: 
<1x4 sparse matrix of type '<class 'numpy.int64'>'
    with 2 stored elements in Compressed Sparse Row format>
In [108]: A.A
Out[108]: array([[1, 0, 2, 0]], dtype=int64)
In [109]: A.T*A           # sparse matrix product, dot
Out[109]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 4 stored elements in Compressed Sparse Column format>
In [110]: _.A
Out[110]: 
array([[1, 0, 2, 0],
       [0, 0, 0, 0],
       [2, 0, 4, 0],
       [0, 0, 0, 0]], dtype=int64)