Special case of sparse matrices multiplication

Question

I'm trying to come up with fast algorithm to find result of operation, where

• L - is symmetric n x n matrix with real numbers.
• A - is sparse n x m matrix, m < n. Each row has one and only one non-zero element, and it's equal to 1. It's also guaranteed that every column has at most two non-zero elements.

I come up with one algorithm, but I feel like there should be something faster than this.

Let's represent every column of A as pair of row numbers with non-zero elements. If a column has only one non-zero element, its row number listed twice. E.g. for the following matrix

Such representation would be

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

Or we can list it as a single array: A = [0, 2, 1, 3, 4, 4]; Now, can be calculated as:

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]]

To calculate we do it one more time:

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]]

The time complexity of such approach is O(mn + mn), and space complexity (to get final result) is O(nn). I'm wondering if it's possible to improve it to O(mm) in terms of space and/or performance?

Answer 1

The second loop combines at most 2m rows of L', so if m is much smaller than n there will be several rows of L' that are never used.

One way to avoid calculating and storing these unused entries is to change your first loop into a function and only calculate the individual elements of L' as they are needed.

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)

This should then have space and complexity O(m*m)

Answer 2

The operation Transpose(A) * L works as follows:

For each column of A we see:

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

The output matrix B = Transpose(A) * L has three rows which are equal to:

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

If we multiply 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)

If you follow through this in a algorithmic way, you should achieve something very similar to what Peter de Rivaz has suggested.

Answer 3

The time complexity of your algorithm is O(n^2), not O(m*n). The rows and columns of L have length n, and the A array has length 2n.

If a[k] is the column where row k of A has a 1, then you can write:

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

and the product, P = A^T*L*A is:

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)

If we turn this around and look at what happens to the elements of L, we see that L[k,l] is added to P[a[k],a[l]], and it's easy to get O(m^2) space complexity using O(n^2) time complexity.

Because a[k] is defined for all k=0..n-1, we know that every element of L must appear somewhere in the product. Because there are O(n^2) distinct elements in L, you can't do better than O(n^2) time complexity.