Особого случай разреженных матриц умножения

Question

Я пытаюсь придумать быстрый алгоритм, чтобы найти результат операций, где

• L - симметричен n x n матрицы с вещественными числами.
• A - редкий n x m матрица, m < n. Каждая строка имеет один и только один ненулевой элемент и равна 1. Также гарантируется, что каждый столбец имеет не более двух ненулевых элементов.

Я придумал один алгоритм, но чувствую, что должно быть что-то быстрее, чем это.

Представим каждый столбец A как пару номеров строк с ненулевыми элементами. Если столбец имеет только один ненулевой элемент, его номер строки указан дважды. Например. для следующей матрицы

Такое представление было бы

column 0: [0, 2]; column 1: [1, 3]; column 2: [4, 4]

Или мы можем перечислить его как единый массив: A = [0, 2, 1, 3, 4, 4]; Теперь может быть вычислена как:

for (i = 0; i < A.length; i += 2): 
    if A[i] != A[i + 1]: 
    # sum of two column vectors, i/2-th column of L' 
    L'[i/2] = L[A[i]] + L[A[i + 1]] 
    else: 
    L'[i/2] = L[A[i]] 

Чтобы рассчитать , мы делаем это еще раз:

for (i = 0; i < A.length; i += 2): 
    if A[i] != A[i + 1]: 
    # sum of two row vectors, i/2-th row of L'' 
    L''[i/2] = L'[A[i]] + L'[A[i + 1]] 
    else: 
    L''[i/2] = L'[A[i]] 

Временная сложность такого подхода является О (м п + т п), и пространство сложность (чтобы получить окончательный результат ) представляет собой О (п п). Мне интересно, можно ли улучшить его до O (m м) с точки зрения пространства и/или производительности?

Answer 1

Вторая петля объединяет не более 2 м строк L ', поэтому, если m много меньше n, то будет несколько строк L', которые никогда не используются.

Один из способов избежать вычисления и сохранения этих неиспользуемых записей - изменить свой первый цикл на функцию и вычислить только отдельные элементы L 'по мере необходимости.

def L'(row,col): 
    i=col*2 
    if A[i] != A[i + 1]: 
    # sum of two column vectors, i/2-th column of L' 
    return L[row][A[i]] + L[row][A[i + 1]] 
    else: 
    return L[row][A[i]] 

for (i = 0; i < A.length; i += 2): 
    if A[i] != A[i + 1]: 
    for (k=0;k<m;k++): 
     L''[i/2][k] = L'(A[i],k) + L'(A[i + 1],k) 
    else: 
    for (k=0;k<m;k++): 
     L''[i/2][k] = L'(A[i],k) 

Это должно затем иметь место и сложность O (м * м)

Answer 2

Операция Transpose(A) * L работает следующим образом:

Для каждого столбца матрицы А мы видим:

column 1 has `1` in row 1 and 3 
column 2 has `1` in row 2 and 4 
column 3 has `1` in row 5 

Выходная матрица B = Transpose(A) * L имеет три ряда, которые равны:

Row(B, 1) = Row(A, 1) + Row(A, 3) 
Row(B, 2) = Row(A, 2) + Row(A, 4) 
Row(B, 3) = Row(A, 5) 

Если мы умножим C = B * A:

Column(C, 1) = Column(B, 1) + Column(B, 3) 
Column(C, 2) = Column(B, 2) + Column(B, 4) 
Column(C, 3) = Column(B, 5) 

Если вы будете следовать через это в алгоритмических образом, вы должны добиться чего-то очень похожее на то, что предложил Питер де Рива.

Answer 3

Сложность времени вашего алгоритма O (n^2), а не O (m * n). Строки и столбцы L имеют длину n, а массив A имеет длину 2n.

Если [к] является столбцом, где строка к А имеет 1, то вы можете написать:

A[k,i] = δ(a[k],i) 

и продукт, Р = А^Т * L * А:

P[i,j] = Σ(k,l) A^T[i,k]*L[k,l]*A[l,j] 
     = Σ(k,l) A[k,i]*L[k,l]*A[l,j] 
     = Σ(k,l) δ(a[k],i)*L[k,l]*δ(a[l],j) 

Если мы повернем это и посмотрим, что происходит с элементами L, мы видим, что L [k, l] добавляется к P [a [k], a [l]], и легко получить сложность O (m^2), используя сложность времени O (n^2).

Поскольку a [k] определено для всех k = 0..n-1, мы знаем, что каждый элемент из L должен появляться где-то в произведении. Поскольку в L есть O (n^2) отдельные элементы, вы не можете сделать лучше, чем O (n^2).