0

I have a large matrix M (M is positive definite) with dimension (10000,10000) and I want to do the following operations:

r = transpose(e_i - e_j) * M * (e_i - e_j)

where e_i is zeros vector(10000 x 1) with all zero entries except at ith entry

I want to do this operation for (e_1,e_2)....(e_1,e_10000)... i.e all pair of i,j belong to {1,10000}.

I have try

-do the calculation directly
-Cholesky factorization

Nevertheless, both methods take too long. They took >=0.5 seconds to finished one calculation and therefore it is not feasible in time.

Is there any methods/library I could try to speed up this process?

Junwei su
  • 138
  • 1
  • 10
  • r = M[i,i] - M[j,i] - M[i,j] + M[j,j] – dmuir Jun 26 '17 at 09:23
  • Possible duplicate of [How to optimize matrix multiplication operation](https://stackoverflow.com/questions/6061921/how-to-optimize-matrix-multiplication-operation) – phuclv Aug 13 '18 at 05:17

2 Answers2

1

You should look into the Strassen Algorithm. It is designed for multiplying large matrices in O(n^2.8). Sample implementation here.

Rivasa
  • 6,510
  • 3
  • 35
  • 64
0

Your expression is transpose(e_i - e_j) * M * (e_i - e_j). Let's take the last part which is M * (e_i - e_j) and multiply it out to get M * e_i - M * e_j.

Now since e_i is zeros except for the ith element M * e_i is just selecting the ith column of M which we can do with no calculation as M[:,i]. Taking the dot product with e_i again this gives us M[i,i].

Multiplying your expression out fully and applying this idea to all 4 terms gives the answer @dmuir gives in the comment above:

r = M[i,i] - M[j,i] - M[i,j] + M[j,j]

This expression should be very fast and does not (in principle) depend on the size of your matrix at all. So you actually don't need any matrix operations at all to compute what you asked for.

Since you said that you want to do this for all i, j pairs you probably want to end up with a matrix R so that R[i, j] is as r given above. You could do that in a loop over i, j or you could vectorise that with something like

d = np.array([np.diag(a)])
R = d + d.T - M - M.T

For the size that you say I would expect this whole operation to complete in around 1 second.

Given the size of the matrices I would probably try to avoid temporaries with

d = np.array([np.diag(a)])
R = d + d.T
R -= M
R -= M.T
Oscar Benjamin
  • 12,649
  • 1
  • 12
  • 14